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ABSTRACT 
Finite  difference  approximations  to   the  one-dimensional  heat 

r       "  ft 

equation  U.  =  U   are  used  to  introduce  explicit  and  implicit  dif- 
ference equations.   The  convergence  and  stability  of  these  equations 
is  discussed  and  these  concepts  are  used  to  show  the  restrictions 
imposed  on  explicit  equations.   The  Implicit  Alternating  Direction 
(IAD)  method  is  introduced  as  a  means  for  solution  of  the  two-di- 
mensional heat  equation.  Although  the  IAD  method  in  its  basic  form 
is  applicable  only  to  parabolic  problems,  it  is  possible  by  slight 
modification  to  apply  the  method  to  elliptic  problems.   Two  examples 
are  used  to  illustrate  the  use  of  the  IAD  method  for  solution  of 
parabolic  and  elliptic  equations  for  a  rectangular  region.  These 
examples  include  a  work  requirement  comparison  with  other  difference 
methods.  A  third  example  is  given  to  show  the  applicability  of  the 
IAD  method  to  non-rectangular  regions,  in  this  case  a  parabolic 
problem  over  a  circular  region.  Results  of  these  examples  show  that 
the  IAD  method  is  a  convenient  and  accurate  method  when  applied  to 
both  parabolic  and  elliptic  partial  differential  equations  and 
suggest  applicability  to  a  wide  range  of  problems. 
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INTRODUCTION 

U   +  U   =  hi.    is  the  parabolic  partial  differential  equation 

governing  the  flow  of  heat  in  two  space  dimensions.  The  constant  k 

appearing  in  this  equation  is  the  diffusivity  of  the  material  through 

which  heat  flows.  We  assume  here  and  subsequently  that  the  time  unit 

is  chosen  such  that  k  =  1.  Numerical  approximations  to  the  solution 

of  this  equation  may  be  obtained  by  solving  an  associated  difference 

equation.  Depending  upon  the  time  level  at  which  the  derivatives 

U   and  0   are  approximated,  the  difference  equation  will  be  one 
xx     yy 

of  two  types:  (1)  explicit,  which  is  simple  to  solve  since  it  involves 
only  one  unknown  quantity  in  terms  of  several  known  quantities,  but 
which  may  also  require  an  enormous  number  of  calculations  due  to 
limitations  on  the  size  of  the  time-step;  and  (2)  implicit,  which 
does  not  restrict  the  size  of  the  time-step,  but  which  involves  five 
unknown  quantities  and  is  most  easily  solved  by  an  iteration  process. 

A  modified  implicit  difference  method,  called  the  Implicit 
Alternating  Direction  (IAD)  method,  was  introduced  by  Peaceman  and 
Rachford  in  1955  and  may  be  used  to  obtain  numerical  approximations 
to  the  solution  of  the  heat  flow  equation  in  a  direct  and  convenient 
manner.  In  addition,  the  IAD  method  may  be  applied,  with  very  slight 
modification  in  form,  to  elliptic  partial  differential  equations. 
The  IAD  method  involves  two  difference  equations:  the  first  in 
which  three  unknown  values  of  U  in  the  x  direction  are  expressed  in 
terms  of  three  known  values  of  U  in  the  y  direction;  and  the  second 
in  which  three  unknown  values  of  U  in  the  y  direction  are  expressed 
in  terms  of  three  known  values  of  U  in  the  x  direction.  Each  of  these 
equations  gives  rise  to  a  system  of  simultaneous  equations.  However, 
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the  coefficient  matrix  of  each  system  is  tri-diagonal ,  and  therefore 
the  solution  is  relatively  easy  to  obtain. 

In  this  paper  we  will  first  consider  the  simple  one-dimensional 
heat  equation  U   =  U..  We  will  derive  the  associated  explicit  and 
implicit  difference  equations  and  then  introduce  the  concepts  of 
convergence  and  stability  to  show  that  a  restriction  on  the  time- 
step  for  the  explicit  equation  is  necessary.  The  explicit  and  im- 
plicit difference  equations  for  the  two-dimensional  heat  flow  equa- 
tion will  be  given  and  the  associated  restrictions  and  solution 
difficulties  will  be  considered. 

We  shall  then  examine  the  difference  equations  which  define  the 

IAD  method,  showing  how  these  equations  are  solved  and  how  they  are 

applied  to  the  elliptic  problem  governed  by  Laplace's  equation 

U   +  U   =  0.  We  will  discuss  some  aspects  of  the  proof  that  the 
xx    yy  r  r 

IAD  method  is  both  convergent  and  stable,  particularly  in  regard  to 
selection  of  an  iteration  parameter  for  use  in  solving  elliptic 
equations. 

Three  examples  are  given.   The  first  example  is  an  application 
of  the  IAD  method  to  the  parabolic  partial  differential  equation  in 
two  space  variables,  and  the  second  is  an  application  to  the  related 
steady  state  problem.   The  region  for  the  problems  in  both  of  these 
examples  is  rectangular.   The  third  example  is  an  application  of  the 
IAD  method  to  the  same  parabolic  equation,  written  in  polar  coordi- 
nates, over  a  circular  region.   In  Examples  one  and  two  a  compar- 
ison is  made  of  the  work  requirement  for  the  IAD  method  and  the  work 
requirement  for  other  difference  methods.   Data  obtained  in  Examples 
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one  and  three  are  used  in  an  attempt  to  obtain  bounds  for  the 
discretization  error. 
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I.     Comparison  of  Explicit  and  Implicit  Difference  Equations  for 
the  Problem  U.    =  U     . 

For  simplicity,   consider  the  parabolic  partial  differential 

equation   (PDE) 

(1)     Ut(x,t)  -  D^Cx.t)     for     a<x<b,       t>0 

with  boundary  and  initial  conditions 

U(a,t)  =  g(t),     U(b,t)  *  h(t),     U(x,0)  =  f(x). 

This  is  the  equation  that  specifies  the  conduction  of  heat  in  one 

space  dimension. 

To  attempt  the  solution  of  this  problem  by  a  finite  difference 

method,  we  establish  a  network  of  grid  points  throughout  the  region 

a^xib,  t  >0.  In  the  x  direction  divide  the  region  into  strips  of 

width  x  =  (b-a)/M,  and  in  the  t  direction  into  strips  of  width 

At  =  T/N,  where  M  and  N  are  arbitrary  positive  integers,  and  T  is 

some  fixed  time  for  which  the  solution  is  desired.  We  may  now  use 

indices  i  and  k  to  denote  the  grid  points  iAx  and  kAt. 

If  we  assume  that  U(x,t)  possesses  a  sufficient  number  of 
partial  derivatives,  then  the  relation  between  U(x,t)  and 
U(x  +ax,  t  +  At)  is  given  by  the  Taylor's  series  expansion 


U(XtAX,t+At)  =    Uix.t)  +  [ax|^  +A-t|l']U6(/0     -f- 


i 
T! 


[ax^  tAtfcfufct)*-  •  •  •  *^L^fc*At3trW> +  R*  > 


where  the  remainder  term  is 

Rn     =     (l/nl)[  A*  |^   +■    *t  ^    ]%(*  +    Ux,   t  +*At) 

for  O^tf^l,  0***1.   This  leads  to  the  following  definition: 

Definition.   If,  in  the  relation  between  U(x,t)  and  U(x  +  Ax, 

t  +  At),  we  neglect  the  remainder  R  ,  then  we  introduce  an  error 

n 
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(called  the  truncation  error)  which  is  of  order  (  Ax  +  At)  ,  denoted 
0(  lx  +  At)  .  Ify  this  we  mean  that  there  exists  a  positive  constant 
C  such  that 

l\\  -  C(  ax  +  At)n  as  both  ax  and  At  go  to  zero. 
Now  if  we  let  v.   ,    be  the  finite  difference  appro ximation  to 

1,K 

U(x,t)  at  the  grid  point  (i  Ax,  kAt),  then  by  using  Taylor's 

expansion  we  obtain  the  following  difference  approximations  for  the 

derivatives  U.    and  U 

t  xx 

vt  =  (vi.k+i "  Ti.k}/  At  +  0(At) 

v»c    =     (vi+l,k-2Ti.k+vi-l,k>^ax>2  +  0<'ix)2- 
Using  these  approximations,  one  finite  difference  approximation  to 
the  PDE  is 

^i.k+l-^.k^**    =     (Vl,k-2Ti,k+vi+lflc)/U3t)2' 
If  we  define  jP  -  (  Ax) /At  and  then  solve  for  v.   .    .   we  get  the 
equation 

<2>  Ti.fcfi=^K.i.k  +  <1-2/)fl>'i.k+(1/^'i+i.k- 

The  boundary  and  initial  conditions  are  approximated  by 
V0,k+1  =  g  r<k+DAt]      •  vMfk+i  -  h[(k+l)At]     ,  and  - 
v .    q  =  f(iAx),  where  i  =  0  and  i  =  M  represent  the  points  x  =  a 
and  X  =  b  respectively,  and  k  =  0  represents  t  =  0.     Equation  (2) 
is  an  explicit  difference  equation  since  we  can  calculate  the  values 
of  v.    .    ^  explicitly  at  time  level  (k  +  l)At  if  we  know  the  values 

of  v.   ,    at  time  kAt.     Thus  by  starting  with  the  boundary  and  initial 
lfK 

conditions,  equation  (2)  can  be  used  to  calculate  the  values  of  v 
at  all  grid  points  for  any  time  t>0. 
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If  we  approximate  the  derivative  U   at  time  level  (&.+  l)At 
rather  than  at  time  level  kAt,  we  obtain  the  difference  approximation 

(vi,k+l  "  vitk)/^t    =    (vi-i,k+l  "  2vi,k+l  +  VLk+lVUx)2. 

2 
Again,  letting  &   =  (  Ax)  /bt   and  rearranging,  we  have 

(3)  vi-i,k+i-(2Y)vi.k+i+vi+i,k+i=  -rvi,k 

Equation  (3)  is  an  implicit  difference  equation  since  it  involves 
three  unknown  quantities  at  time  level  (k+l)At  in  terras  of  one 
known  quantity  at  time  kAt.  The  boundary  and  initial  conditions 
have  the  same  approximations  as  for  the  explicit  equation.   Thus 
we  have  a  more  complicated  scheme  in  this  case  since  at  any  time 
level  (k+l)A>t,  the  implicit  equation  (3)  is  written  once  for  each 
point  l^iiM-1,  resulting  in  the  following  system  of  M-l  simul- 
taneous equations  in  the  M-l  unknowns  v.  .  .  , 

( i  +  2/)*KKt(  -  yf  vZj Ktl  '     =  ^  +  yf  ^DAt] 

Now  let  &±  =  o±  =  -l/tf       and  b±  =   (H-2/tf  ),    (i  =  1,2,..., M-l), 
and  call  the  known  terms  on  the  right  hand  sides  of  these  equations 
d-.,  d2,   ....,  cL.  _ .     Then  after  dropping  the  subscripts  k+1,  the 
system  has  the  following  form 

Vl  +  ClV2  =     ^L 

a2Vl  +  b2V2  +  C2V3  =  d2 


aiVi-l  +  biVi  +  Vi+1         ■  di 


"M-2M-3  +  ^-2^-2  +  ^-2^-1  =  ^-2 
^-1^-2  +  Vl^-l      =  ^-1 
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The  matrix  of  coefficients  for  this  system  is  tri-diagonal ,  and 
the  use  of  Gaussian  elimination  results  in  the  following  algorithm 
for  solution. 

Let     /$,=  b,   >      r,=  d,^3( 

then      \ZM.,  =  Vi      )       V^rrlQ^fi)^'    )     (i-M-Z.M-V  '  -,l)  . 

Thus  we  have  our  choice  of  two  finite  difference  schemes  with 
which  we  can  attempt  to  find  an  approximate  solution  to  the  PEE  (l). 
We  shall  see  that  one  of  the  conditions  required  to  make  the 
approximation  a  useful  one  is  a  restriction  on  the  size  of  the  time- 
step  At  when  using  the  explicit  equation  (2).  It  will  be  shown  in 
Example  1  that  the  use  of  the  implicit  equation  (3)  requires  less 
work  than  the  use  of  (2),  even  with  the  added  complexity  of  solving 
systems  of  equations  since  there  is  no  restriction  on  the  size  of  At 
in  the  implicit  method. 

Whether  or  not  the  solution  v.  .  of  the  difference  equation  is 
a  good  approximation  to  the  solution  U(x,t)  of  the  PDE  can  be  answered 
by  considering  the  questions  of  convergence  and  stability.  We  begin 
with  the  following  definitions: 

Definition.  The  difference  w.  .  ■  U(i  ax,  k  At)  -  v.  .  is  called 
the  discretization  error  at  the  point  (iAx,  kAt). 
Definition.  We  say  that  the  solution  v.  .  of  the  difference  equation 

"" — — — — — ■  2.  ,K 

converges  to  the  solution  U(x,t)  of  the  PDE  at  the  point  (i  Ax,  kAt) 

if  lim    w.  .  =  0.  If  v.  .  converges  to  U(x,t)  at  every  grid  point 
AtjAtoO  X»K  1»K 

(iix,  k^t),  then  we  say  that  the  difference  procedure  isnconvergent. 
When  we  speak  of  convergence,  we  are  speaking  about  the  difference 
between  the  exact  solution  to  the  PDE  and  the  exact  solution  to  the 
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difference  equation.   Since  we  will  generally  attempt  to  solve  the 

difference  equation  using  a  digital  computer,  we  cannot  expect  to 

obtain  the  exact  solution.   Instead  we  obtain  a  numerical  solution 

which  differs  from  the  exact  solution  due  to  round-off  error.  When 

ihis  difference  remains  small,  we  say  that  the  corresponding  procedure 

is  stable. 

Before  giving  a  formal  definition  of  stability,  we  examine 

the  question  of  round-off  error  in  more  detail.     Suppose  v(x,t)  is 

the  exact  solution  to  the  difference  equation  and  suppose  we  replace 

v(x   ,   t  )  by  v(x   ,t   )  +  e  at  the  grid  point  (x   ,   t   ).      This  may  be 

considered  to  be  the  effect  of  round-off  error  or  possibly  a  mistake 

in  the  calculations.     We  call  e  the  error  at  the  point   (x  ,   t  ).      If 

o   o 

we  then  continue  the  calculation  using  the  value  v(x  ,  t  )  +  e,  with- 
out introducing  new  errors,  we  obtain  a  solution  v*(x,t)  at  the 
point  (x,t).  We  call  the  difference  v*(x,t)  -  v(x,t)  the  departure 

of  the  solution  caused  by  the  error  e  at  (x  ,  t  ).   If  errors  are 

J  o   o 

committed  at  more  than  one  point,  we  speak  of  the  cumulative  depar- 
ture due  to  these  errors.   For  linear  problems  we  can  use  the 
principle  of  superposition  to  conclude  that  the  departure  due  to  two 
or  more  errors  is  the  sum  of  the  departures  due  to  the  individual 
errors.  We  can  now  give  the  following  definition: 
Definition.  A  difference  procedure  is  stable  if  the  cumulative 
departure,  due  to  errors  in  the  solution,  tends  to  zero  as  the 
maximum  of  these  errors  tends  to  zero  and  does  not  "grow  faster"  than 
some  positive  power  of  1/ax  as  ax  tends  to  zero.  We  say  that  the 
cumulative  departure  does  not  "grow  faster"  than  some  positive  power 
of  l/ax  to  exclude  exponential  growth  in  1/ax. 
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The  general  method  used  to  show  convergence  is  to  find  a  bound 
for  the  discretization  error  w(x,t).  and  then  show  that  this  bound 
depends  on  ax  and  At  and  tends  to  zero  as  ax,  At->0.  Stability 
is  most  often  proved  by  a  method  due  to  Von  Neumann.   This  method 
involves  expanding  the  error  in  a  Fourier  series  and  showing  by  means 
of  the  coefficients  of  this  series  that  this  error  does  not  grow  as 
the  solution  progresses.   For  the  two  difference  schemes  which  we 
have  introduced,  we  will  obtain  the  conditions  for  convergence  and 
stability  by  a  method  given  in  Forsythe  and  Wasow  L  3 J  . 

We  assume,  for  convenience,  that  the  associated  boundary  condi- 
tions for  the  PDE  (l)  are  g(t)  =  0  and  h(t)  =  0.  Without  loss  of 
generality  we  can  assume  that  a  »  0  and  b  =  tt  since  we  can  modify 
this  interval  by  means  of  a  linear  transformation.  'Then  the  exact 

solution  to  (1)  is 

4?       2 
(k)     U(x,t)  =   2   a  e"r  t  sin  rx, 

where  the  initial  function  f(x)  has  been  expanded  in  the  convergent 

Fourier  sine  series 

ex?  fir 

f(x)  -        /*  .   a  sin  rx,  where  a  =  ^r    f (x)  sin  rx  dx. 
Now  consider  the  following  difference  equation 

«)  (\,k+1  -\,k>/*t  =  *'(*i+i.k+l-2vi.k+l  +  'i-l.k+l^4*)2 
+  0<r)   (v1+1>k-2v1>k  ♦  v^^/C  Ax)2. 

where  o~  is  a  parameter.  Note  that  for  cr-  -  ot  equation  (5) 
reduces  to  the  explicit  difference  equation  (2);  and  for  O"   =1,  it 
reduces  to  the  implicit  equation  (3).  We  can  now  examine  the  stabil- 
ity and  convergence  of  equations  (2)  and  (3). 
Theorem.  When  cr  <  1/2,  equation  (5)  represents  a  convergent  and 
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p 
stable  difference  scheme  if  V?=  ( Ax)  /a  t  ^  2(l-2cr).  When 

c7_  — 1/2,  equation  (5)  represents  a  convergent  and  stable  difference 

scheme  for  all  values  of  vtf  . 

Proof.   Suppose  that  v(x,t)  is  a  solution  to  (5)  which  satisfies  the 

conditions:  v(x,0)  =  f(x),  0  <  x  <  IT   ;  v(0,t)  =  v(tt  ,t)  =0,0  <  t  ^T. 

Here,  for  convenience  of  notation,  we  use  (x,t)  to  mean  the  grid  point 

(i  a.x,  kAt).   Assume  that  f(x)  can  be  written  in  the  Fourier  series 

oa 

a  sin  rx,  where 
r 


*ie   now  seek  a  solution  of  the  form 


v(x,t)  =    y   ,  ar  7  r(t)  sin  rx. 
If  y   (0)  =  1,  the  solution  satisfies  the  initial  condition;  the 
boundary  conditions  are  automatically  satisfied.   If  we  substitute 
y   (t)  sin  rx  into  (5).  we  get 

With  the  initial  condition  y   (0)  =  1,  this  difference  equation  has 
the  so^ 


r 
■lution    t/,^)     -     <P** 


H--      ±3 


3-  VAX 


where         (D    _    I  -  &0' *)*'*    *#* 


+  4  &  5i^L  ^|X 

and  where  t/At  takes  only  positive  integer  values.   So  we  have  as 

a  solution  <%> 

Z^nt/At 
arY^r  sin  rx, 

V-l 

CK> 

if  this  series  converges.   Since    /       |"-nH  ^     ,  a  sufficient 
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condition  for  convergence  (in  fact,  for  uniform  convergence)  is 
\Y   \  £     1  for  all  r.   But  this  means 

The  right  hand  inequality  is  always  satisfied.   The  left  hand 
inequality  is  satisfied  when  2(1  -  2<r)^  v0  .  Note  that  this  is 
always  the  case  when  0~^l/2,  so  that  \y  \    -      1  for  any  ]<?  . 

To  show  convergence  it  is  necessary  to  show  that  the  series 

(6)  tends  to  that  given  by  (4)  as  At,Ax->  0. 

CD  ?         2     2 

In  the  expression  for  y     we  replace  sin  r  ax  by  r  (ax)  (h 

.  V  2. 

+  0(Ax)  to  get 


#  = 


I  +  A t  [cr  r z  4  4<T  0O*)x] 


We  now  examine  the  Lim     Y  .   If  we  consider  this  as  an 

iterated  limit,  we  have  two  cases.   These  are: 

Casel-    kSS^o^    =  ttlftL- hf^tte — J 

•      I  At-90l  |4-tr\r2At    J         )  ~ 

Case  2.        Lim       Lim      Y  .     As  it  stands,  Lim  (P     '  is 

an  indeterminate  form  so  we  write      l_ikvi     ytr        =     LJyvi  CT 

Using  L' Hospital's  rule  we  have 

t  Lm     ±!$  =  I  Urn    *U  <Q.  =  t.[+G-l°-)0(**)-*1]. 
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So     Urn  Uw   %        =    Um    tf  =  e 

Thus  we  see  that  regardless  of  the  way  in  which&x  and^t  tend 
to  zero,  Lim         (D  =     e"  so  that  every  term  of  (6) 

tends  to  the  corresponding  term  of  (b). 

Since  the  convergence  of  (6)  is  uniform,   it  is  possible  to 
consider  the  limit,   as  Ax, A  t  tend  to  zero,   in  a  termwise  fashion. 
Therefore  if  2(1  -  2<r)f  ^  ,  then  lim  v(x,t)  ■  U(x,t),   and  the 
difference  equation   (5)  represents  n   convergent  scheme. 

To  show  stability  it  is  more  convenient  to  represent  the 

values  of  v(x,t)  at  each  grid  point  by  means  of  a  finite  Fourier 

series  rather  than  the  infinite  series   (6).      This  is  based  on  the 

theory  of  trigonometric  interpolation  and  uses  the  orthogonality 

relations 

M-l 

M-l  fl-i 

2L>  cos     nrax  =        /iL>     sin     nrAx  =  M/2,  0<n<M, 

where  M  =  T/a  x  is  an  integer.      Using  these  relations,  we  can 
write                          M., 

f(x)  =         x    ,  A     sin  nx, 


n-/ 


n 


where  M.( 


A       =(2/M)     /   ,         f(rAx)   sin  nrAx. 


)7L 

This  relation  will  be  valid  for  the  points  x  =  rAx,(r  =  1,2,. . .  fM-y. 

We  use  a  finite  series   since  round-off  error  cannot  be  regarded  as 

a  smooth  function  independent  of    Ax. 

Now,   let  M  =   t/ax  be  the  number  of  grid  points  on  a  line  for 
constant  t  and  consider  an  arbitrary  line  of  errors  e  ,(p=l,2,. . . ,M-l) 
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at  t  =  0.     We  then  consider  these  errors  as  some  initial  function 

and  represent  them  by 

M-l 


e       =  ^__,  A     sin  sr  ax 

S   /  ™  n  r 

(s  =  1,2,..., M-l) 


where  M_, 

A       =    (2/M)      ^L     e     sinprAx. 
r  p-i        P 

(r  =  1,2,..., M-l) 

Then  by  the  same  argument  that  resulted  in  the  infinite  series   (6), 

we  show  that  ^  . 

(7)  v*(sAx,nAt)     =     2^   Arji       sin  sr^x, 

\r=/ 

is  the  solution  of  (5)  at  the  points  of  the  grid  for  the  initial 

values  e    ,  or  v*  represents  the  departure  due  to  e     for  the  line  t  =  0. 

We  now  determine  a  bound  for  v*.      Substituting  the  expression 

for  A     in  (7)  we  have  . .  , 

r  M-l 

(8)  v*(s  Ax,nAt)  =  2L    Sp.s,nep 

Where  M-l  t/Lt 

gp  s  n     =   (2/M)         /    ,     YV      sin  srAx  sin  prax 

(Here,  our  interchange  of  summations  is  valid  since  these  are 

finite  series).      If  2(1  -  2  cr)  £  kD  ,  then     \%\  4    1     and    |g  |  ±  2. 

Let     |e    If  a    f°r  r  =  lt2,...,M-l  where     0  represents  the  maximum 


r 
of  the  errors.      Then  from   (8)  we  have 

v*(s6x,nAt)     £     2M&        =       2  "HT  5    /ax 

which  means  that  the  cumulative  departure  due  to  errors  is  of  order 

o(ax)~  which  by  our  definition  means  that  the  process  is  stable. 

Therefore,  we  have  proved  that  for  ^£2(1  -  2 CT) ,   equation  (5) 

represents  a  convergent  and  stable  finite  difference  scheme.   It 
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follows  then  that  for  <J~  >   1/2,  (5)  represents  a  convergent  and 
stable  scheme  for  any  positive  value  of  \0   •   For  cT  =  0,  (5) 

reduces  to  equation  (2)  and  the  condition  for  convergence  and 

2 

stability  is  j^^2,  or  At  £  (Ax)  /2.   As  previously  stated,  the 

explicit  method  is  restrictive  in  the  size  of  the  time-step.   For 
CJ~=  1.  (5)  reduces  to  the  implicit  equation  (3)  and  is  unrestricted 
as  to  the  size  of  the  time-step  allowed. 

The  method  used  to  prove  stability  in  the  previous  theorem  is 
essentially  the  method  of  Von  Neumann.  This  method  can  be  extended 
to  problems  involving  two  or  more  space  variables. 

Note  that  for  the  difference  equations  examined,  the  conditions 
for  stability  and  convergence  are  identical.  Lax  and  Richtmyer  j_5J 
have  developed  a  theory  for  a  large  class  of  linear  partial 
differential  equations  with  constant  coefficients  which,  for  slight- 
ly different  definitions  of  convergence  and  stability,  states  that 
stability  of  a  difference  equation  implies  convergence  of  the  equa- 
tion.  Briefly,  this  theory  assumes  that  the  initial  functions  of 
the  linear  differential  problem  belong  to  some  Banach  space  B,  of 
functions  of  x.   Then  associated  with  every  initial  function  in  B 
is  a  solution  to  the  difference  problem  (this  solution  being  a 
function  of  x  alone  when  t  is  fixed)  as  well  as  a  so-called  genuine 
solution  to  the  PDE.   Both  of  these  functions  belong  to  B.   Under 
these  conditions  a  finite  difference  procedure  is  called  convergent 
if  the  solution  to  the  associated  difference  equation  converges  to 
a  genuine  solution  for  all  initial  functions  in  B.   This  is  a  more 
restrictive  definition  than  the  one  we  have  used  since  we  have 
established  convergence  for  initial  functions  in  some  incomplete 
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subset  of  the  space  B.  A  difference  procedure  is  called  stable  in 
this  sense  if  the  solution  v(x,t),  corresponding  to  some  initial 
function  f(x),  satisfies 

||v(x,t)||  £     M(t)  ||   f(x)||   ,  0  i  t  i  T,  for  all  f(x)  in  B, 
where  M(t)  is  independent  of  Ax.  This  definition  is  also  more 
restrictive  than  ours  since  we  allow  a  bound  which  depends  on  A  x. 
There  is  also  a  difference  in  the  concept  of  distance  between  two 
functions.  We  have  used  the  norm 

llv(x.t)  -  U(x,t)l|   =  sup  |v(x,t)  -  U(x,t)|   , 
whereas  Lax  uses  the  mean  square  norm 

If  the  definitions  of  convergence  and  stability  are  taken  in  the 
above  sense ,  then  convergence  and  stability  are  proven  to  be  equiv- 
alent concepts  in  the  Equivalence  Theorem  of  Lax.   It  can  be  shown 
that  explicit  difference  equation  (2),  with  the  restrictions  on  At, 
actually  does  satisfy  the  Lax-Richtmyer  theory. 
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II.   The  Explicit  and  Implicit  Difference  Equations  for  the  Problem 
in  Two  Space  Variables. 

For  the  problem  U   +  U   =  U.  ,  we  use  the  same  difference 
xx    yy    t' 

approximations  for  U   and  U.  and  a  similar  approximation  for  U 

xx     x>  yy 

(the  approximations  for  U   and  U   taken  at  time  level  k^t).  We 
^  xx     yy 

change  notation  slightly,  letting  the  subscripts  i  and  j  refer  to 
increments  in  the  x  and  y  directions,  and  the  superscript  k  refer 
to  increments  in  time.   The  explicit  difference  equation  is 
(9)   ^j  =     VhJ    f  p  L^-'o  +  H  +  V'  +  v/0JH-  4^jJ 

Again  taking  the  approximation  for  U   and  U   at  time  level 
(k  +  1)  At  results  in  the  implicit  difference  equation 

do)   v-  ♦  <; , „-  t „«-,  -  (4,  fM™  -  -  f  tf 

As  in  the  case  of  the  difference  equations  in  one  space  vari- 
able, we  find  that  the  explicit  procedure  (9)  is  both  stable  and 
convergent  for  certain  values  of  &    and  the  implicit  jorocecfure  (10) 

is  stable  and  convergent  for  all  values  of  ^ .   (Here  we  have  taken 

2  2 

J£  =  (  Ax)  /^t  =  (  Ay)  I  At,)     The  restriction  on  f   for  the 

explicit  equation  is  f  '  ^»   This  result  is  proven  in  Douglas 
|_6 J  ,  in  which  the  difference  equations  are  treated  as  matrix 
operators  so  that  the  result  is  valid  for  any  number  of  space  dimen- 
sions. 

2 
For  the  explicit  equation,  V?£  k,   means  that  At  £.    (ax)  /k, 

so  that  the  allowable  time- step  is  half  the  size  of  that  allowed 

for  the  explicit  equation  in  one  space  variable.   This  may  mean  that 

an  enormous  number  of  calculations  are  required  to  obtain  a  solution 

for  a  particular  time  T. 
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liquation  (10)  again  results  in  a  system  of  simultaneous 
equations  which  can  be  solved  by  Gaussian  elimination.   However, 
such  a  method  for  solution  is  impractical  as  there  are  now  five 
unknowns  per  equation  and  the  number  of  required  computations  is 
more  than  twice  the  number  required  for  the  tri-diagonal  system. 
The  solution  of  this  equation  is  possible  by  iteration  techniques , 
but  we  shall  demonstrate  the  solution  of  the  PDE  in  two  space  vari- 
ables by  means  of  a  modified  implicit  method  called  the  Implicit 
Alternating  Direction  (IAD)  method. 
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III.   The  IAD  Method 

We  now  introduce  the  IAD  method  of  Peaceman  and  Rachford  [7  J  • 

If  only  one  of  the  second  partial  derivatives  for  the  PDE,  U   +  U 

^  xx    yy 

=  U+  f  say  U  ,  is  replaced  by  a  difference  approximation  at  time 
(k  +  l)At,  and  U   is  replaced  by  an  approximation  at  time  level 
kAt,  then  we  again  have  a  difference  equation  with  only  three  unknowns, 
and  this  equation  may  be  solved  by  the  use  of  the  previously  stated 
algorithm.   The  difference  equation  obtained  is  said  to  be  implicit 
in  the  x  direction.  We  then  repeat  the  process  for  a  second  time- 
step  of  equal  size,  reversing  the  level  at  which  the  second  deriva- 
tives are  approximated,  to  obtain  a  difference  equation  implicit  in 
the  y  direction.  When  writing  these  equation,  we  use  time  levels 
2k  +  1  and  2k  +  2  to  emphasize  that  the  time-steps  must  be  of  equal 
size.   The  equations  are: 
Implicit  in  the  x  direction 

C11}    *V*.  Tax? o^ft 


Implicit  in  the  y  direction 
(12) 


XK*i  iM|  y2Ktl  2KH         ,2-K+l  ZKtZ.  tT.Ktl        Kil. 


AVi  (ax)1-  (a^)1 

2  2 

We   choose    ax  =   Ay,   and  let      V?  =   (Ax)   /a  t  =  (  Ay)   /At.      Then, 

rearranging,  we  obtain  a  form  more  convenient  for  calculation: 
(11a)      -VL_h}+Z(f+\Hl}     -Vcil/j         =    ^-i+i(f'Kj    +JijJ*i 

(12a)      -^j-i+Uf+'M.j      -^,       =     ^j**(f-'M/     tU^0J 

Each  of  these  equations  leads  to  M  -  1  equation  in  M  -  1  un- 
knowns at  each  half  time-step,  At/2.   The  systems  each  have  a  tri- 
diagonal  coefficient  matrix. 
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In  addition  to  the  time-dependent  problem  U       +  U       =  U .  ,  the 
IAD  method  may  also  be  used  in  the  solution  of  steady  state 
problems;   for  example,   Laplace's  equation,  U       +  U       =0.      The  pro- 
cedure is  essentially  one  of  iteration  with  each  stage  of  iteration 
being  regarded  as  a  time-step  in  a  time -dependent  problem.     A  start- 
ing value  is  used  which  approximates  the  boundary  conditions  and 
this  serves  as  an  initial  condition.     Equations   (11a)  and   (12a)  are 
used  as  alternate  stages  of  iteration  with    ^  serving  as  an  iteration 
parameter.     We  will  show  by  means  of  examples  how  the  IAD  method  is 
used  in  the  solution  of  both  time-dependent  and  steady  state  prob- 
lems. 

We  now  state  the  following  theorem  without  proof. 
Theorem.      The  IAD  method,   as  applied  to  the  partial  differential 
equations    AU(x,y,t)  =  U, (x,y,t)  and    AU(x,y)  =  0,  having  either 
Dirichlet  or  Von  Neumann  boundary  conditions,   is  both  convergent 
and  stable  for  any  region. 

The  proof  of  this  theorem  can  be  found  in  Forsythe  and  Wasow 

[3]    ,  pp  272  and  411,   and  in  Birkhoff  and  Varga    [&]    .     We  now 

state  some  of  the  results  of  the  proof  given  in  [3J      for  the  IAD 

method  as  applied  to  the  steady  state  problem. 

Consider  the  difference  approximations  for  U       and  U       used  in 

xx     yy 

equations  (11)  and  (12).  If  we  regard  the  solution  of  the  difference 
equation  as  a  column  vector  V(x,y),  then  the  symmetric  difference 
expression  used  to  approximate  U   can  be  regarded  as  a  matrix  A 
operating  on  V.   Thus, 

AhV(xfy)  =  -V(x-ax.y)  +  2V(x,y)  -  V(x+Ax,y). 
Similarly,  for  U  ,  the  analogous  matrix  is  A  ,  and 
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AVV(x,y)  =  -V(x.y-Ay)   +  2V(x,y)   -  V(x,y+^y), 

2 
For  the  case  of  Laplace's  equation  in  a  square  region  with  n 

u  0  0 

interior  grid  points,  A  will  be  an  n  x  n  matrix  which  may  be 

partitioned  into  an  n  x  n  block  diagonal  matrix  with  each  diagonal 

v 
block  being  an  n  x  n  tri -diagonal  matrix.  A  will  be  similar  to 

such  a  matrix.   If  we  use  the  subscript  k  to  denote  the  iteration 

number,  the  IAD  method  is  then  defined  by  the  equations 

(13)    VK.,4  --  V*.  -£(M*«.  *  AVi) 

w    vK    =  W-£(aX*  +  AVJ 


where  the  parameter  ^-  tt-  may  vary  from  iteration  to  iteration. 
In  effect,  we  are  solving  the  elliptic  (steady  state)  problem  as  the 
asymptotic  solution  of  the  time -dependent  problem;  thus  the  At 
appearing  in  the  expression  for  s#  .   In  the  actual  solution  process 
for  the  elliptic  problem,  however,  bt   will  not  enter  into  the  method 
for  determining  optimal  values  of  ^fi . 

We  now  define  the  error  E,  as  the  difference  V,  -  U,  where  U 
is  the  solution  of  the  PDE  at  each  grid  point  and  may  thus  be  re- 
garded as  a  vector.   Then,  after  eliminating  V„  1  /?  between  (13) 
and  (1^) ,  the  effect  of  one  complete  iteration  is 

where  the  matrix  operator  P  is 

pit)  =  (i-  ^t' ri-  ^  a")(i + ^  Ahr '(r- 

and  I  is  the  identity  matrix. 

Now,  as  shown  by  Forsythe  and  Wasow,  a  sufficient  condition  for 
convergence  is  that  the  spectral  radius  of  P(j^)  be  less  than  one, 
where  the  spectral  radius  is  the  maximum  of  the  absolute  values 
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of  the  eigenvalues  of  P.  This  is  shown  to  be  true  for  any  positive 
constant  value  of  y?. 

If  the  region  of  interest  for  the  PDE  is  a  rectangle,  spec- 
ific statements  can  be  made  about  the  values  of  the  iteration  para- 
meter ^?.  which  lead  to  the  most  rapid  convergence.  As  shown  in 
Birkhoff  and  Varga  L8J  »  however,  this  is  not  possible  in  non- 
rectangular  regions.   Suppose,  for  convenience,  that  the  boundary  of 
the  region  is  a  square  of  side  length  TT  .   Then  A  and  A  commute 
and,  by  a  theorem  due  to  Frobenius,  both  matrices  share  the  eigen- 
vectors 

sin  px  sin  qy,  (p,q  =  1,2 M-l), 

with  eigenvalues 

2  2 

4  sin  pdx/2,  ^  sin  qAy/2,  (p,q  =  1,2,. . .  ,M-l). 

The  eigenvalues  of  P(P  )  are  then 

a*  %::s$}8;is#}  ><>*'» "-o 

From  this  expression  it  can  be  seen  that  the  eigenvalues  of  P(^) 
are  less  than  one  in  absolute  value  for  any  positive  value  of  "f . 
Peaceman  and  Rachford  j_7j  use  the  Von  Neumann  method  of  stability 
analysis  to  derive  the  same  expression  as  an  amplification  factor 
for  error.  In  Example  2  we  will  show  how  this  expression  is  used 
to  determine  the  optimal  values  of  the  iteration  parameter  u? . 

We  will  show  in  Examples  1  and  2  how  the  IAD  method  can  be 
used  to  obtain  an  approximate  solution  to  the  time-dependent  and 
the  steady  state  problems  for  a  rectangular  region.   In  Example 
3  the  method  will  be  applied  to  time-dependent  problem  over  a 
section  of  a  circular  region. 
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rjcample  1.      Solution  of  a  Time-Dependent  Problem  Using  the   IAD 
Method. 

Problem:      To   find  the  temperature  distribution  for  various 
times   in  a  square  plate  with  edge-length  one. 

The  temperature  distribution  in  a  square  plate  of  edge-length 
one  (See  Figure  24)  is  governed  by  the  following  parabolic  partial 
differential  equation: 

^(x.y.t)   +  UyyCx.y.t)     =   UtCX;y,tJ   , 
with  boundary  conditions 

U(l,y,t)     =     U(x,l,t)     =     0 
and  initial   condition 

U(x,y,0)  =  1 


U^  o 


Ux(0,y,t)  =  U  (x.O.t)  =  0 


>  x 


Figure  24  Temperatures  in  a  Square  Plate 
To  solve  this  problem  using  the  IAD  method  we  use  the  differ- 


ence equations 


,iMI         ,ZMl 


IK 


(12a)       -\/C)yi     ♦  JLCf+O'tj     ■  ^jt,      =     /fro  j     +ity-»Kj     W«.;J 


^   _      (ax)1     _      <£  ^ 
where  y    -      —gj: 


&  t 
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Equation  (lla)  is  said  to  be  Implicit  in  the  x  direction  and 
equation  (12a)  is  said  to  be  illicit  in  the  y  direction. 

To  apply  these  equations,  we  divide  the  region  shown  in  Figure 
2k   into  N  increments  of  length  l/N  in  both  the  x  and  y  directions. 

Here  N  is  an  arbitrary  positive  integer.   The  resulting  grid  is  shown 

2k+l 
in  Figure  25.   Thus  v.  .  refers  to  the  difference  approximation  to 

*•  J 

the  solution  U(x,y,t)  at  the  grid  point  (iAx,  j  Ay)  at  time  (2k+l)  At  , 
where  Ax  =  Ay  =  l/N. 

if 


CO;MA^) 


(0;j6,d) 


(rJAX,  Nfcj) 


OJAXjjAij) 


Ci*x,o) 


^AX,o) 


Figure  25  Grid  for  the  Region  Shown  in  Figure  2k 
Equations  (lla)  and  (12a)  can  be  used  for  all  interior  points 
of  the  grid  and,  of  course,  the  values  of  U  are  known  on  the  boundar- 
ies x  =  1  and  y  =  1.  However,  at  the  boundaries  x  =  0  and  y  =  0,  the 
boundary  conditions  are  U  (0,y,t)  =  0  and  U  (x,0,t)  -  0,  so  we  must 
find  an  approximation  involving  U  for  U   at  x  =  0.  and  an  approx- 
imation involving  U  for  U   at  y  =  0.   Since  the  procedure  is  the 

y    yy 

same  at  both  boundaries,  the  details  are  given  for  the  boundary  x  =  0 
only. 

We  assume,  for  a  fictitious  point,  a  distance  ax  outside  the 
boundary,  that  U(x-  &x,y,t)  =  U(x+  Ax,y,t).   This  is  intuitively 
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acceptable  since  U  =0  and  the  region  under  consideration  can  be 
thought  of  as  one  quarter  of  a  square  plate  having  edge-length 
two  and  being  centered  at  the  origin.   For  this  larger  plate  (with 
boundaries  held  at  temperature  zero)  there  would  be  complete  sym- 
metry and  hence  no  heat  flow  across  the  x  and  y  axes.  We  now  write 
the  Taylor  expansions  for  U(x-  Ax,y,t)  and  U(x+  Ax,y,t)  in  terms  of 
U(x,y,t),  where  (x,ytt)  represents  the  boundary  point 
[o,j  Ay,(2k+1)  A  t]  .   Thus  we  have 
U(x-AX,y,t)  =  U(x,y,t)  -  AxU^x.y.t)  +  ^U^x.y.t) 

"  ^Uxxix^t>  + 

and  t 

U(x+  Ax.y.t)  =  U(x,y,t)  +  AxU(x,y,t)  +  (£*)  U  (x,y,t) 

X  n   |    XX 

A* 


+   ^IlU*^t)  +  — 


Now,  adding  these  equations,  and  neglecting  terms  involving  U 

AAaA 

and  all  higher  derivatives,   we  get 

U(x-  Ax,y,t)  +  U(x+Ax,y,t)  =  2U(x,y,t)   +  (  Ax)2U     (x,y,t)   +  0(  A  x)' 

But  U(x-  Ax,y,t)  =  U(x+ax,y,t),   so  we  can  write 

(16)  Uxx(x,y,t)  =  ^[uU+AX.y.t)   -  U(x,y,t)j      +  0(  Ax)2. 

(17)  Uyy(x,y,t)  =  (A*[u(x,y+ Ay.t)  -  U(x,y,t)]     +  0(  Ay)2. 
For  notation  in  the  computer  program  IADMT,  we  use  TST(l,J)  to 

denote  the  values  of  our  approximation  to   U(x,y,t)   at  the  end  of  a 
first  half  time-step    (using  equation  11a)  and  TEND(I,J)   to  denote 
the   approximation  at  the  end  of  a  second  half  time-step    (using  equa- 
tion 12a).      Thus   TKND(I,J)   represents   the  numerical   solution  at  the 
end  of  a  complete   time-step.      In  the   equations   that   follow,   we 
abbreviate   TST(I.J)   by  T*(I,J)   and   TEND(I,J)   by  T(I,J). 
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All  oomputer  programs  were  run  on  the  IBM  3&0  computer  using  the 
FORTRAN  IV  language.   Since  FORTRAN  IV  does  not  allow  zero  subscripts 
we  must  increase  our  indices  by  one  and  thus  take  N  =  11  to  obtain 
ten  increments  of  length  Ax  =  Ay  =  0.1.   Then,  T(1,J),  (J=l,2,. . . ,11) 

corresponds  to  v2^1  ,  ( j=0,l, . . . ,10),  and  T(I,l),  (1=1,2 11) 

corresponds  to  v.  +  ,  (i=0,l,. . . ,10) ,  at  time  (2k+l)  Zkt. 

Using  this  notation,  and  the  approximations  (16)  and  (1?),  the 
difference  equations  at  the  boundaries  x  =  0  and  y  =  0  are 

T*(1.J)  -  T(l.J)  u     2T*(2.J)  -  2T*(1.J)   T(l. J-l)-2T(l.J)+T(l. J+l) 


At/2       "      (Ax)2        + (  Ay)2 


or 


(lib)  2(*+l)T*(l,J)-2T*(2,J)  =  T(l,J-l)+2(f  -l)T(l,J)+T(l,J+l) 
for  (J=2,3 N-l)  ' 

2(^+l)T*(l,l)-2T*(2,l)  =  2T(l,2)+2(/>-l)T(l,l) 
'  for  J  =  1, 


and 


or 


T(I.l)-T»(I.l)  _  T*(I-l.l)-2T*(I.l)+T*(I+l.l)   2T(I.2)-2T(I.l). 
4t/2      =       (ax)2  Uy)2 


(12b)  2(y*+l)T(I,l)-2T(I,2)  =  T*(I-l,l)+2(f  -l)T*(I,l)+T*(l+l,l) 
for  (1=2,3 N-l) 

2(vo+l)T(lfl)  =2T(1,2)  =  2T*(2,1)  +  2(y?-l)  T*(l,l) 
'    for  I  =  1  ' 

Now  using  equations  (11a),  (lib),  (12a),  and  (12b)  we  can  write 

out  the  systems  of  simultaneous  equations  which  are  to  be  solved. 

Theses  are  I 

Implicit  in  the  x  direction 

Z(fH)T*(W)-ZT*U)J)  =  D(1) 

-t*o,t;+^ot*<*j?-t*(3,,t; 

-Tm(*lJT)+l(f*i)r*(*i,Jl-T*(*hJ)~     ~     ~~-  $-2) 
-T(i*'l,T)  +  Z(f*i)T*(i4'tjT)  =  D(m-|) 
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where  D(I)  e  T(I.J-l)  +  2(v* -1)T(I,J)  +  T(I,J+l) 
(1=1,2 N-i)    (J=2,3 N-l), 

and   D(I)  =  2T(I,2)  +  2(]?-l)  T(I,l) 
(1=1,2 N-l)      J=l 

Implicit  in  the  y  direction 
Z(f+i)T(I,i)-lT(l,Z)  =     D(l) 

-T(I,l)+^ftOT(r^)-TCl,3)  =  D(2) 

-T(r/~l)+Z(f*Or(I,J)-T(r,Jtl)  =  D(J) 

-TCl.H-V  +  Uf+nTd^-D-TCt^-t)  =  D(N-2) 

-Td^-l^iCffOrCl^-i)  =  D(N-l), 

where  D(J)  =  T*(I-1,J)  +  2(f -1)T*(I,J)  +  T*(I+1,J) 
(J  =  1,2,...,N-1)     (I  -  2,3 N-l) 

and    D(J)  =  2T*(2,J)  +  2(^-l)T*(l,J) 

(J  =  1,2 N-l)     1=1. 

If  we  call  the  coefficients  of  the  main  diagonal  terms  B(I) 

and  the  coefficients  of  the  lower  diagonal  and  upper  diagonal  terms 

A(I)  and  C(I)  respectively  »  then 

A(I)  =  -1,  (I  =  2,3 N-l)  C(I)  =  -1,  (I  =  1,2 N-2) 

B(I)  =  2(^+1),  (I  =  1,2,..., N-l).   In  addition,  let  F  =  2(^-1), 

The  algorithm  for  solving  these  systems  is  then 

/SO)'  BCD  j  *f0)  =   D0)/fi(i) 

/5(i)=  B(I)- A(i)C(i-0/fl(i'O  ,     ^Ci)^i)-Ki)^(i-\)]/f5(L) 

(I  =  2,3 N-l) 

r*(N-i)=rcM-/)  ,   T*ti)=rw-c(JiT*(in)/flu) 

77N-i)  =  r04-i)  ,     raj  =  TS(z)  -CCDTCt+O/fia) 

(I  =  N-2.N-3 1) 
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The  flow  charts  for  the  solution  program  (fcADMT)  and  the  algo- 
rithm subroutine   (TDIAG)  are  given  in  Figures  1  and  2.      Twelve  time- 
steps,     At  =  0.001,   0.001,   0.001,  0.002,   0.005,  0.01,  0.01,  0.02, 
0.05,   0.05,   0.05,   and  0.05,  were  used  consecutively  to  obtain  solu- 
tions at  times  0.001,   0.002,   0.003,   0.005,   0.01,   0.02,   0.03,   0.05, 
0.10,  0.15,  0.20,   and  0.25, 

The  analytic  solution  for  the  problem  was  obtained  by  the 
separation  of  variables  technique,  and  is 

U(x,y,t)     =   (l6/Tr2)ZZS%^T)€  Cosa^NjXcotfn-iJW 

A  program  was  written  to  obtain  the  analytic  solution  at  the  grid 
points   (i  ax,   j  Ay),    (i,j  =  1,2.,. . .  ,N-l),  for  the  times  0.005,  0.01, 
0.05,   and  0.250.      This  solution  is  called  ANSOL(l.J),  and  includes 
all  terms  in  the  double  summation  for  n,m  =  1  to  20.      The  values 
of  ANS0L(I,J)  were  read  into  the  program  IADMT  for  comparison  with 
the  values  of  the  numerical  solution,   TEND(I,J).      The  difference, 
W(I,J)  =  ANS0L(I,J)  -  TEND(I,J),  and  the  percent  difference, 
Z(I,J)  =  100  W(I,J)/ANS0L(I,J),  were  calculated  for  all  points  with 
the  exception  of  those  points  on  the  boundaries  where  the  temper- 
ature is  held  at  zero. 

Figure  3  shows  a  print-out  of  the  numerical  solution  and  analy- 
tic solution  for  time  0.25.      Figure  k  shows  the  difference  and  per- 
cent difference  for  the  same  time. 

Numerical  solutions  were  obtained  for  several  values  of    At 
and    ax,   all  at  time  0.005.      Table  1  lists  the  values  of  the  numer- 
ical solution     /(x,y,t)     and  the  percent  difference.    Z(I,J).   at  the 

point  (0.6,0.6).     The  analytic  solution  for  this  point  is  VCx^ft)    -. 
=  0.9999. 
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Table  1.      Numerical  Solution  and  Percent  Difference  at  Point 
(0.6.0.6)   for  Different  Values  of    U  and    Ax. 


At 

AX 

(*x)2 

At  +  (  Ax)2 

v(x.y.t) 

z(i,J) 

0.0005 

0.025 

0.00062 

0.00112 

0.9996 

0.033 

0.001 

0.025 

0.00062 

0.00162 

0.9995 

0.044 

0.0025 

0.025 

0.00062 

0.00312 

0.9991 

0.088 

0.005 

0.025 

0.00062 

0.00562 

0.9981 

0.180 

0.0005 

0.05 

0.0025 

0.0030 

0.9994 

0.055 

0.001 

0.05 

0.0025 

0.0035 

0.9994 

0.055 

0.0025 

0.05 

0.0025 

0.0050 

0.9990 

0.099 

0.005 

0.05 

0.0025 

0.0075 

0.9982 

0.170 

0.0005 

0.1 

0.01 

0.0105 

0.9975 

0.240 

0.001 

0.1 

0.01 

0.0110 

0.9975 

0.240 

0.0025 

0.1 

0.01 

0.0125 

0.9973 

0.260 

0.005 

0.1 

0.01 

0.0150 

0.9965 

0.340 

All  solution  programs  were  run  on  the  IEM  3&0  computer.      This 
computer  uses  seven  digits  to  the  right  of  the  decimal  point  in  a 
decimal  number.      Therefore,  we  can  write  as  a  bound  for  the  round- 
off error  C  x  10      ,  where  C  is  a  positive  constant.      (See  McCracken 
and  Dorn     [10 J    . )     C  is  determined  by  the  number  of  arithmetic 
operations  performed  and  by  how  the  round-off  accumulates  during 
these  operations.      It  is  therefore  possible  for  C  to  become  quite 
large.      However,   the  satisfactory  results  obtained  in  all  solutions 
lead  us  to  believe  that  random  cancellation  of  round-off  has  occurred, 
thus  keeping  C  small.      Furthermore,   since  the  discretization  error 
depends  on    At  and   (  ax)    ,   and  the  smallest  of  these  values  was 
6. 2  x  10     ,  we  shall  neglect  the  round-off  error  and  examine  only 
the  discretization  error. 

From  Douglas     £9]    *   the  discretization  error  is  bounded  by  the 

2 
quantity  C,(At)  +  CjAx)    ,  where  C,    and  C?  are  positive  constants 

and  depend  on  the  least  upper  bounds  of  the  derivatives   U. .    and 

U  .      We   can  write  this  bound  as  K    -At  +  (Ax)     (  ,   where  K  =  maxCn  ,CL, 

xxxx  "-  J  12 

Thus  we   say  that  the  discretization  error  is  of  0  [_^t  +  (Ax)     J  . 
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For  convenience,  we  have  used  At  x  l(r  and  (  Ax)  x  l(r  when 
working  with  the  data  from  Table  1.   This  data  was  used  to  plot 
percent  error  versus  L  At  +  (ax)  J   x  10  ;  the  result  is  shown  in 
Figure  5.   The  line  shown  in  this  figure  was  obtained  by  a  least 
squares  fit.  Although  the  data  appears  to  fit  a  straight  line,  we 
can  not  take  this  as  verification  that  the  discretization  error  is, 
in  fact,  of  0  [_At  +  (  ^x)  J  .   To  obtain  verification  we  would 
need  to  know,  or  at  least  have  some  estimate  of,  the  derivatives 
U^.,.  and  0...  Knowing  the  analytic  solution,  it  would  be  possible 
to  calculate  these  derivatives;  this  was  not  done  since  the  objective 
of  the  paper  is  to  illustrate  the  use  of  the  IAD  method  rather  than 
to  verify  predicted  results. 

We  now  compare  the  work  required  using  the  IAD  method  against 
that  required  when  using  the  explicit  method.  We  make  this  com- 
parison for  the  solution  at  time  0.25.  Although  the  IAD  method  is 
stable  for  any  size  time-step,  we  have  seen  that  the  discretization 
error  depends  on  At.   Therefore,  we  avoid  using  very  large  time- 
steps.   For  this  problem,  0,25  was  the  final  value  of  t,  and  solutions 
were  obtained  for  11  smaller  values  of  t  by  using  values  of  At  rang- 
ing from  0.001  to  0.05.   To  find  the  number  of  arithmetic  operations 
required  for  the  IAD  method  in  this  problem,  we  take  as  our  start- 
ing point  the  solution  of  equation  (11a),  which  gives  the  temper- 
atures implicit  in  the  x  direction.  As  we  have  seen,  this  equation 
results  in  N-l  =  10  equations  in  N-l  unknowns.   The  algorithm  used 
to  solve  these  equations  involves  three  expressions,  each  of  which 
required  three  arithmetic  operations  (including  addition  and  subtrac- 
tion).  These  expressions  were  formed  a  total  of  N-2  times,  so  that 
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9(N-2)   arithmetic  operations  are  required  for  the   solution  of  each 
system  of  N-l  equations.      To  obtain  the  temperature  at  every  grid 
point  at  the  end  of  a  half  time-step,  equation   (11a)  must  be  used  N-l 
times.      Thus   for  a  half  time-step,   9(N-2)*(N-l)  operations   are  requir- 
ed.     Now,   the   solution  of  equation   (12a),  which  gives  the  tempera- 
tures implicit  in  the  y  direction,   requires  the  same  number  of  oper- 
ations for  a  half  time-step.      Thus  the  solution  for  a  complete  time- 
step  requires  2x9(N-2)»(N-l)  operations.      As   stated,   a  solution  at 
time  0.25  was  obtained  as  the  twelfth  complete  time-step,   so  that  the 
total  number  of  arithmetic  operations  required  using  the  IAD  method 
was 

12x2x9 (N-2) -(N-l)     =     1?,W 

2 

Use  of  the  explicit  equation   (9)   requires  that    At  t=.   (ax)   fh  s 

0.0025.      This  means  that  the  solution  must  be  stepped  through  100 
time-steps  to  reach  time  0.25.      The  solution  of  equation  (9)   in- 
volves 7  arithmetic  operations,   and  this  equation  must  be  solved  at 
each  of  the   (N-l)2  grid  points.     A  total  of  100x7 (N-l)2  =  70,000 
operations   are  required   for  a  solution  at  time  0.25.      Thus  we   see 
that  the   IAD  method  requires  two-thirds  less  work  than  the  explicit 
method  for  this  problem.      This  problem  involved  100  grid  points.     As 
the  number  of  grid  points  increases,   it  becomes  even  more  advantageous 
to  use  the  IAD  method  rather  than  the  explicit  method.      Note  also 
that  the  maximum  time-step  used  for  this  problem  was  0.05,   which  is 
20  times  larger  than  the  maximum  allowable  time-step  for  the  explicit 
method. 
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Example  2.   Solution  of  a  Steady  State  Problem  Using  the  IAD  Method. 

Problem:  To  find  the  steady  state  temperature  distribution  at 
the  face  of  an  infinite  prism  with  edge-length  one. 

The  temperature  distribution  for  this  infinite  prism  (See 
Figure  26)  is  governed  by  the  elliptic  partial  differential  equation 

Uxx(x»y)  +  Uyy(x,y)  =  °* 
with  boundary  conditions 

U(0,y)  =  U(x,0)  =  0, 

U(x,l)  =  1 
Ux(l,y)  =  -hU(l,y). 
The  last  boundary  condition  means  that  the  boundary  x  =  1  is  a 
partially-conducting  boundary.   For  convenience  we  take  the  value  of 
the  constant  h  to  be  one. 

In  solving  this  problem  by  the  IAD  method,  we  treat  it  as  the 
limiting  case  of  the  time-dependent  problem  U   +  U   =  U. .   This 
allows  us  to  use  the  same  difference  equations,  (11a)  and  (12a), 
that  we   used  in  Example  1, 


U=0 


»  UX  =  -^U 


U  =  o 


>  x 


Figure  26.   Temperatures  at  the  Face  of  an  Infinite  Prism. 
We  divide  the  region  shown  in  Figure  26  in  exactly  the  same  way 
in  which  we  divided  the  region  for  Example  1.   Thus  the  resulting 
grid  is  identical  to  the  grid  shown  in  Figure  25. 
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Tdn,   IAD  method,  as  applied  to  this  steady  state  problem,  can 
be  considered  as  a  process  of  iteration.   In  such  a  process  we  use 
a  starting  approximation  to  the  solution,  and  then  attempt  to  improve 
this  approximation,   We  call  each  attempt  at  improvement  an  iteration. 
Thus  we  regard  the  use  of  both  equations  (11a)  and  (12a)  as  one  com- 
plete iteration,  and  we  regard  the  quantity  y?  appearing  in  these 
equations  as  an  iteration  parameter.  We  now  wish  to  determine  the 
optimum  values  of  y?  ,  i.e.  ,  the  value  or  values  of  V3  which  will  give 
us  the  most  improvement  in  our  starting  approximation  after  the  small- 
est number  of  iterations. 

rtecall  from  the  discussion  of  the  convergence  of  the  IAD  method 

that,  for  a  rectangular  region,  it  is  possible  to  predict  optimum  values 

of  \Q  .   In  that  discussion,  the  effect  of  one  complete  iteration  was 

represented  by  E.  =  P(  ^?,  )E,  ,  , 

where  P(^,  )  is  a  matrix  operator  and  K  is  the  difference  between 

the  true  solution,  U(x,y),  and  the  numerical  solution,  v  (x,y), 

obtained  after  the  kth  iteration.  We  call  E,  the  error  at  the  kth 

k 

iteration.   We  saw  previously  that  we  were  assured  convergence  (for 
an  infinite  number  of  iterations)  if  the  eigenvalues  of  the  operator 
P  were  all  less  than  one  in  absolute  value.   These  eigenvalues  are 
given  by, 

(l5a)  0°-4Slo1-p£?)(y-4S.H^%fr) 

(The  difference  between  this  expression  and  (15)  is  due  to  the  fact 
that  the  region  for  this  problem  has  side-length  one). 

From  (15a)  we  see  that  each  of  the  N-l  eigenvalues  of  P  can  be 
reduced  to  zero  on  some  one  of  N-l  iterations  if  we  use  values  of 
P    piven  by 
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(18)  ^p  =  4  sin  pTT/2N,  (p  =  1,2 N-]). 

Thus,  theoretically  we  should  be  able  to  reduce  the  error  to  zero 
after  N-l  =  10  iterations.   Of  course,  complete  error  reduction  in 
10  iterations  could  happen  only  when  using  an  exact  numerical  pro- 
cedure.  Due  to  round-off  error  we  do  not  have  an  exact  procedure, 
so  we  can  only  hope  that  error  will  be  reduced  to  a  very  small,  and, 
therefore,  acceptable  level.   Moreover,  since  in  the  IAD  method,  one 
complete  iteration  involves  the  solution  of  both  equations  (Ha) 
and  (12a),  we  require  20  intermediate  iterations  for  maximum  error 
reduction. 

Using  equation  (18),  the  values  of  &   and  the  associated  values 
of  At  (from  At  =  (ax)  /s£  )  were  calculated  and  are  listed  in 
Table  2. 

Table  2.   Calculated  values  of  the  iteration  parameter  \0 

At 

0.40612 
0.04587 
0.0170? 
0.00916 
0.00593 
0.00432 
0.00344 

0.00293 
0.00264 

0.00252 

Since  this  is  a  steady  state  problem,  At  has  no  real  meaning.   We 
use  it  here  only  to  emphasize  the  similarity  in  the  time-dependent 
and  steady  state  problems  when  using  the  IAD  method.   As  in  the  time- 
dependent  problem,  where  we  obtained  solutions  in  order  of  increas- 
in|   At,  we  also  use  the  values  of  \0  given  in  Table  2  in  order  of 
increasing  At,  i.e.,  we  use  the  larger  values  of  ^f>   first.   Each 


&3 


-E- 

p. 

1 

0.02462 

2 

0.21799 

3 

0.58579 

4 

1.09202 

5 

1.68713 

6 

2.31287 

7 

2.90^98 

8 

3.41421 

9 

3.78201 

10 

3.97538 

value  of  ^  must  be  used  over  a  complete  iteration  so  that  any 
particular  y2  is  used  for  two  intermediate  iterations. 

It  is  possible  to  reduce  the  number  of  iterations  by  grouping 
values  of  s£>  .   From  Table  2  it  can  be  seen  that  for  larger  p  the 
difference  between  values  of  V?  is  small.   Therefore,  an  average 
value  of  y?  may  be  used  for  certain  groups  of  ^P   .   The  eigenvalues 
of  P  corresponding  to  the  values  of  p  falling  within  that  group, 
while  not  zero ,  are  sufficiently  small  so  that  effective  error 
reduction  still  occurs.   For  reasonably  small  N  (as  in  this  case), 
it  is  possible  to  calculate  all  values  of  y?  and  set  up  the  groups 
by  inspection.   A  procedure  for  setting  up  groups  when  N  is  large 
is  given  in  [7]  .  For  this  problem  we  group  y?  as  in  Table  3. 

Table  3.   Grouped  Values  of  V?  . 


_P_ 


avgAt  f=   (Ax)2/A  t 


8,9,10  0.00270  3.7037 

5,6,7  0.00455  2.1978 

3 ,4  0.01310  0.7633 

2  O.04587  0.2179 

1  0.40612  0.0246 

Thus,   by  grouping  values  of  f>  ,  we  are  able  to  reduce  the  num- 
ber of  iterations  required  by  one-half  with  the  expectation  that  the 
results  obtained  should  be  nearly  as  good  as  those  obtained  when 
using  all  calculated  values  of  y?  .      Douglas  and  Peaceman  LllJ      give 
examples  in  which  simple  values  of    s^  are  used.      They  start  with  a 
value  between  4  and  10,   and  then  halve  this  value  repeatedly  to 
obtain  additional  values  of  ^P  .      We   shall  obtain  solutions   to  this 
problem  using  all  three   approaches. 

For  a  starting  approximation,   corresponding  to  the  initial 
condition  in  the  unsteady  state  problem,   we  use 
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(19)  U(x,y)  =  xy/(l  -  y  +  x). 

This  expression  satisfies  the  boundary  conditions  at  all  but  the 
conducting  boundary. 

As  in  Example  1,  we  can  use  equations  (lla)  and  (12a)  at  all 
interior  grid  points  and  at  the  boundaries  where  the  values  of 
U(x,y)  are  known.  At  the  conducting  boundary,  we  may  use  the  same 
approximations  for  U,  and  U   but  must  find  a  new  approximation  for 
U  .   Again  we  assume  a  fictitious  point,  a  distance  ax  outside  the 
boundary,  at  which  the  temperature  is  U(x+Ax,y),  and  as  before  we 
write  the  Tay3.or  expansions  for  U(x+Ax,y)  and  U(x- Ax,y)  in  terms 
of  temperatures  U(x,y)  on  the  boundary,  where  (x,y)  represents  the 
point  [N  Ax.jAyJ  .   Combining  these  expressions,  we  have,  as  before 

(20)  U(x-  Ax,y)  +  U(x+ Ax,y)  =  2U(x,y)  +  (  Ax)2U  (x,y)  +  0(ax)2. 
However,  we  no  longer  have  the  s\jmmetry  which  allows  us  to  say  that 
U(x-  Ax,y)  =  U(x+  Ax,y).  We  must  use  the  boundary  condition  to 
evaluate  U(x+  Ax,y).   The  method  here  is  one  due  to  Schmidt  as 
given  in  Carslaw  and  Jaeger  L12J  .  We  use  a  symmetric  difference 
for  the  normal  derivative  at  the  boundary  to  get 

Ux(x,y)  =  l/fcAx)[u(x+Ax,y)  -  U(x-  Ax,y)]   +  0(Ax)2 
oolving  for  U(x+Ax,y),  and  using  the  boundary  condition 

Ux(x,y)  =  -U(x,y) 
we  have 

U(x+Ax,y)  =  U(x-Ax,y)  -  2AxU(x,y)  +  0(  Ax)2. 
Substituting  in  (20),  we  have,  finally, 

(21)  U^Cx.y)  =  2/(  Ax)2  [  U(x-  Ax,y)  -  (1+  Ax)0(x,y)]   +  0(  Ax)2. 
In  the  computer  program  IAEM  we  use  TST(I.J)  to  denote  the 

approximation  to  U(x,y)  at  the  end  of  each  intermediate  iteration 
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and  TEND(I.J)  to  denote  the  approximation  at  the  end  of  a  complete 

iteration.   In  the  equations  that  follow  we  abbreviate  TST(I,J)  by 

T*(I,J)  and  T3ND(I,J)  by  T(I,J).   Using  this  notation  and  the 

approximation  (21),  the  difference  equations  at  the  conducting 

boundary  are : 

Implicit  in  the  x  direction 

T»(N.J)  -  T(N.J)  =  2T*(N-1.J)  -  2(1+  Ax)T*(N.J) 

(  Ax) 

+  T(N.J-l)  -  2T(N.J)  +  T(N.J+l) 

(Ay)2 
or 

(lie)  -2T*(N-1,J)  +  [^(k)  +  2(1+ ax)]  T*(N,J) 

=  T(N,J-l)  +  Lf(k)   -  z]   T(N,J)  +  T(N,J+l), 

where 

f(k)     =  (Ax)2/Atk  =  Uy)2/*^ 

and  k  denotes  the  iteration  number. 

Implicit  in  the  y  direction  : 

T(N.J)  -  T*(M.J)     2T*(N-1.J)  -  2(1+  Ax)T*(N,J) 

(Ax) 

+  T(N.J-l)  -  2T(N.J)  +  T(N.J+1). 


(Ay)2 


or 


(12c)  -T(N,J-1)  +  [f(k)   +  2]   T(N,J)  -  T(N,J+1) 

=  2T*(N-1,J)  +  [f(k)   -  2(1+ Ax)]  T*(N,J). 
We  empahsize  here  that  the  value  ^(k)  must  be  used  for  a  complete 
iteration  before  introducing  the  value  ^(k+l). 

Now,  using  equations  (11a),  (12a),  (lie),  and  (12c)  we  write 
out  the  systems  of  equations  to  be  solved: 
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Implicit  in  the  x  direction 
0  +  yW]-ra,J-)  -  T*(3,J)  =     d(2) 

-T*(2^)+J^^^K)]T*(3;I)-T*^J)  =     D(3) 

-TVx,J)  +&+  W]TV'^)  -  T*fN,J)  =     D(N-l) 

-T*(M-lJJ)+-[f^^O^K)jT*Cfs/,j)      =     D(n) 

(J=2,3 N-l) 

where       D(I)  =  T(I,J-l)  +    [^(k)  -  2]  T(I,J)  +  T(I,J+l) 

(1=2,3 N)  (J=2,3 N-l) 

As  in  Example  1,  the  coefficients  are  A(l),  B(I)  and  C(I),  with 

B(I)  =  2  +  f(k),    (1=2,3 N-l),  B(N)  =f(k)   +  2(1+ Ax) 

A(I)  =  -1,  (I=3A..., N-l),  A(N)  =  -2 

C(I)  =  -1,  (I  =  2, 3,..., N-l) 
Implicit  in  the  y  direction 
lZ+f(K)]T(I,l)  -  7(1,1)  =  D(2) 

-TU.iULi+fMtfCi.y-Tii^    __    __  =    D(3) 

-TCij-^  +  litfOxilTdjj)  -  T(r,Jv<)  =    d(j) 

-T(l,N-3)  +  [z->fO\)lT(l,M-Z)~  r(I,M-l)  -     D(N-2) 

-T(I,M-2)  +  \X+f(*)]T(.Z,H~\)         =     D(N-l) 

(1=2,3 N) 

where     D(J)  =  T*(I-1,J)  +    ty(k)  -  2]  T*(I,J)  +  T*(I+1,J) 

(1=2, 3,..., N-l)  (J-2,3 N-2) 

D(N-l)  =  T*(I-1,N-1)  +    fy>(k)  -  2]  T*(I,N-1)  +  T*(I+1,N-1) 
+  T(IfN) 
(1=2, 3,..., N-l) 
and     D(J)  =  2T*(N-1,J)  +    [f(k)   -  2(1+ax)]    T*(N,J) 
I=N  ,    (J=2,3 N-2) 
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D(N-l)  =  2T*(N-l,N-l)  +  [f(k)   -  2(1+ ax)]  T*(N,N-1)  +  T(N,N). 
We  note  that  D(N-l)  contains  the  term  T(I,N)  which  has  been  trans- 
ferred from  the  left-hand  side  since  it  is  a  known  quantity,  by- 
virtue  of  the  boundary  condition  U(x,l)  =  1. 

The  coefficients  are  written 

B(I)  =  2  +  f(k),   (1=2,3 N-l) 

\(I)  =  -1,    (1=3,4 N-l) 

C(I)  =  -1,    (1=2,3 N-2)   . 

We  use  the  algorithm  given  in  Example  1  to  solve  these  two 
systems  of  equations.  The  flow  charts  for  the  solution  program 
IAUM  and  the  algorithm  subroutine  TDIAG  are  given  in  Figures  6  and  2. 

The  analytic  solution  to  this  problem  was  obtained  from  Churchill 
l_13j   and  is 


oo 


(22)  U(x,y)  =        *Z  ^^^J sin  «  x 


where       A   =  1  -  cos  0( 
n     n 

1  +  cos2o( 


n 


and  (X-,  ,  0(^  ,  . . . .  are  the  roots  of  tan  o(  =  -  £(. 
^  program  was  written  to  obtain  the  analytic  solution  at  the  grid 
points  (iAx,jAy),  (i,  j=2,3,. . .  tN) ,  and  this  solution  is  called 
V(I,J).   To  obtain  V(I,J),  the  series  in  (22)  was  summed  until  a 
term  was  reached  for  which  the  absolute  value  of  that  term,  divided 
by  the  sum  of  previous  terms,  was  less  than  10  .   The  summation  was 
terminated  upon  reaching  such  a  term.   The  values  of  V(I,J)  were 
read  into  the  program  IAEM  to  allow  comparison  with  the  numerical 
solution,  TEND(I.J).   IADM  calculates  the  difference, 

//(I, J)  =  V(I,J)  -  TEND(I.J)  and  the  percent  difference, 
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Z(I,J)  =  100  W(I,J)/V(I,J)  for  all  grid  points  except  at  the 
known  boundaries,  x  =  y  =  0  and  y  =  1. 

The  first  solution  was  obtained  using  the  five  grouped  values 
of  y?  listed  in  Table  3.  Figures  7  through  16  give  the  results  for 
all  five  iterations.   In  these  figures  as  well  as  in  those which 
follow,  the  odd  numbered  figures  show  the  numerical  and  analytic 
solutions  for  a  particular  iteration  while  the  even  numbered  figures 
show  the  difference  and  percent  difference.   Results  of  all  five 
iterations  are  shown,  so  that  the  effect  of  each  value  of  y?  and 
the  process  of  error  reduction  can  be  seen.  We  note  that  error 
reduction  takes  place  diagonally  from  the  larger  grid  points  to  the 
smaller.   This  is  because  the  larger  values  of  y2  which  are  used 
first  correspond  to  the  larger  eigenvalues  of  the  operator  P,  and 
these  eigenvalues  in  turn  correspond  to  larger  values  of  x  and  y. 

The  second  solution  was  obtained  using  the  ten  values  of  Y? 
listed  in  Table  2.  Results  of  the  tenth  iteration  are  shown  in 
Figures  17  and  18.   Comparison  of  Figures  l6  and  18  shows  that  the 
results  of  using  the  grouped  values  of  Y?  in  five  iterations  are 
essentially  the  same  as  the  results  of  using  all  calculated  values 
of  V2  in  ten  iterations.   In  both  cases  errors  range  from  approxi- 
mately 0.05  percent  to  2.5  percent.  We  therefore  conclude  that  by 
grouping  values  of  y2  we  can  reduce  the  number  of  required  iterations 
significantly. 

The  third  solution  was  obtained  using  the  five  simple  values 
of  \0  :  4,2,1,0.5,  and  0.25.   Results  are  shown  in  Figures  19  and  20. 
It  can  be  seen  that  error  reduction  was  less  than  in  the  solution 
obtained  by  using  five  grouped  values  of  yO  .   However,  these  results 
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are  considered  to   be   satisfactory  when  the   resulting  simplification 
of  the  problem  is  taken  into  account. 

We  use   the  results  of  the  second  solution   (Figure  16)   to  make 
a   comparison  of  the  work  requirement  for  the   IAD  method  and  the  work 
requirement  for  another  iteration  method.      In  the   IAD  method,   the 
equation  U(x,y)  =  xy/(l  -  y  +  x)  was  used  to  obtain  starting  values 
at  each  grid  point.      Comparison  of  these  values  with  the   analytic 
solution  shows  a  maximum  difference   (maximum  initial   error)  of 
approximately  0.10.      Comparison  of  the  numerical  solution  for  the 
last  iteration  and  the  analytic  solution  shows  a  maximum  difference 
(maximum  final  error)  of  0.00*1-1.      Thus  the  reduction  of  error  is  on 
the  order  of  4x10     .      In  Example  1  we  found  that  a  complete  time- 
step  using  the   IaD  method  required  2x9(N-2)* (N-l)   arithmetic  oper- 
ations  (including  addition  and   subtraction).      For  this  problem  a 
complete  time-step   corresponds   to   a  complete   iteration,   and  since 
five   complete  iterations  were  used ,   the   total  number  of  operations 
required   for  an  error  reduction  of  4x10    "  was   5*2x9 (N-2)» (N-l) 
=  8100.      Peaceman  and  Rachford    [7 J    give  the  work  requirement  for 
the  extrapolated  Liebmann  iteration  method  of  solution  to  Laplace's 

equation  as  5v(N-l)     for  an  error  reduction  of  10       .      Thus,   for 

-2 
an  error  reduction  of  10      ,    the  number  of  arithmetic  operations   by 

the  extrapolated  Liebmann  method  is  5x2 (N-l)     =  10.000.      Thus  the   IAD 

method  requires   fewer  operations   for  a  larger  error  reduction  than 

the  extrapolated  Liebmann  method,    considered  to  be  one  of  the  more 

efficient  methods  of  iteration. 
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Example  3.   Solution  of  a  Time -Dependent  Problem  in  a  Non-rectang- 
ular Region  Using  the  IAD  Method. 

Problem:   To  find  the  temperature  with  time  for  a  quarter  section 
of  the  face  of  an  infinite  cylinder  with  radius  one. 

The  temperature  at  the  face  of  this  infinite  cylinder  (See 
Figure  27)  is  governed  by  the  following  parabolic  partial  differ- 
ential equation 

(23)  u   +— u  + -L ,  uQA  =  -i-u. 
rr   r    r   ra  00    k  t 

with  boundary  and  initial  conditions 

u(r,0,t)  =  u(r,  TT/2,t)  =  1,  u(l,0,t)  =  1  and  u(r,9,0)  =  0 

As  before,  we  choose  the  time  unit  such  that  the  diffusivity  k  is  one. 
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Figure  27.      Temperatures  at  the  Face  of  an  Infinite  Cylinder. 
We  divide  the  region  shown  in  Figure  27  into  N  increments  of 
length    Ar  =  l/N  in  the  r  direction,   and  into  N  angles  of  magnitude 
£0  =  Tr/2N  in  the  0  direction.      The  resulting  grid  is  shown  in  Figure 
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Figure  28.      Grid  for  the  Region  Shown  in  Figure  27. 
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Let  v(r,e,t)  be  the  difference  approximation  to  u(r,6,t).  We 

then  use  the  following  difference  approximations  for  the  partial 

derivatives  u , ,  u  ,  u   ,  and  iUQ: 
u   r   rr      oo 

v,f ^—^ +  0C^) 

Now  consider  the  two  terms  u       +(l/r)u   ,  of  (23).      They  are  approxi- 
mated by 

y(^ANT,9,t)-  */(V-A/-,a,-t)      +      y(V- &r, e;-t) - z y/(v, e> t) + \/(W&^ e,t ) 

2-^ANT  Cav^ 

<a^t-  f   UC^J 

The  complete  difference  equation  approximating  (23)  is 
(24)   V(^6/fcfAt)-  ^,0,*) 

-[Ofr*)  +0(^)L4  0(^)1]    =  0. 


Now  consider  (24)  as  an  operator  L  (v)  =0.   If  L  operates 
on  u(r,0tt),  we  have 


iW 
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where   the  tilde  indicates  that  the  derivatives  are  to  be  taken  at 
certain  points  in  the  region  0^r<l,   0<6    <TT/2,   0£t  £T. 

We  call  z(r,0,t)  =  v(r,0,t)   -  u(r,9,t)   the  discretization  error. 
Then  since  L     is   a  linear  operator,   we  have 

Lr(v)  -Lr(u)  =  Lr(v  -  u), 
or 

Lr(z)    =  -fuit  r  IffifafXm*  +  &$**«>*  +  0(M)\  0<*rAcW 

Call  the  right  hand  side  of  this  equation  w(r,0,t).   Thus  we  see 
that  z(r,0,t)  is  a  solution  to  the  difference  equation 

Lr(z)  =  w(r,9,t) 
Since  (24)  represents  a  convergent  difference  scheme,  we  have 
that 

l.u.b  (z|  ^  l.u.b  (  w|  , 
or 

(25)   |z|  6.    CxAt  +  C2(  Ar)2  +  C3(  Ad)2, 

where  CU  ,  C0  and  Co  depend  on  the  least  upper  bounds  of  u. . ,  u 

i       c.  j       *  rir  ut,       rrrr 

and  ufiftflfi  respectively.      In  addition,   C0  and  C~  depend  on  r.      There- 
fore, we  should  expect  a  discretization  error  of  0   [_At  +  (  Ar)   +(&&)]. 

We  now  digress  somewhat  to  consider  the  IAD  method  as  applied 
to  the  related  elliptic  PDE  (u  +(l/r)u  +(l/r  )  vlqq  =  0)  for  this 
non-rectangular  region.  If  U(r,8)  represents  a  vector  solution  to 
this  equation,   then  the  difference  analog  to    (24)  is 

~£7HAh  +  AVJ      U(r,6)     =     0 
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h       v 

where  A  and  A  are  square  matrices  such  that 

AhU(r,6)  =  -   (I-   f£)  U(r-Ar.e)  +  2U(r,0)  -  (1  +  $)U(r+ar,9) 

AVu(r,0)  = "     ~^i^[?irt9~*e) " 2U(r,e)  +  U(r»6+A6)]    • 

For  the  simple  case  of  N  =  4  (a  net  with  9  interior  points)  we 


find  that 
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h  v 

For  this  case,  although  A     and  A     are  non-singular,  neither  is 

symmetric,   and  a  calculation  shows  that  A  A     ^  A  A  .      Thus,   for 

h  v 

this  non-rectangular  region,  A     and  A     do  not  share  a  common  eigen- 
vector basis,   and  it  would  therefore  not  be  possible  by  the  method 
used  in  Example  2  to  predict  values  of  the  iteration  parameter  for 
optimum  convergence  rate.      It  would  be  necessary  to  use  simple 
values   for  this  iteration  parameter  (See  Example  2).      This  is  con- 
sidered to  be  the  most  serious  disadvantage  of  the   IAD  method. 
However,   for  the  parabolic  problem  under  consideration,  we  shall 
see  that  entirely  satisfactory  results  are  possible  even  though 
the  region  is  non-rectangular. 
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We  now  return  to  equation  (24).      Call  V.     .     the  approximation 

1 » J 

to  u(r,e,t)  at  the  grid  point     (iAr,   j  a6,  kAt)    .      Let    f>  =  (*r)  /At 
and    (X  =  A  6;   also  note  that  r  =  iAr.      Then  equation   (24)  yields  the 
"working  equations"  analogous  to  equations    (lla)   and   (12a).      These 
are: 

Implicit  in  the  r  direction 


Implicit  in  the  0  direction 

(27)     -  jfel  ^j-,   +  Zfo+fjl'c.j     ~  TFk--  %t, 


2JC-+I  ,    .       ,     ,£K-H         ,  ■     s       IK+ 


Since  the  coefficients  in  these  equations  depend  on  i,  this 
will  add  some  complexity  to  the  computer  program.   However,  since 
the  temperature  is  given  at  all  boundaries,  we  need  find  no  special 
approximations  at  any  boundary,  as  was  necessary  in  Examples  1  and  2. 
We  use  the  same  notation  used  in  the  previous  two  examples  when 
writing  out  the  two  systems  of  simultaneous  equations  which  result 
from  Equations  (26)  and  (27).   These  are  : 

Implicit  in  the  r  direction 
2(^>+l)T*(2,J)  -  (1+1/2 )T* (3, J)  =  D(2) 

-(l-l/4)T*(2,J)+2(f  +l)T*(3,J)-(l+l/4)T*(4,J)  =  D(3) 

-(1-1/2(1-1) )T*(I-l,J)+2(^+l)T*(I,J)-(l+l/2(I-l))T*(I+l,J)=D(I) 

-l(l-l/2(N-3))T*(N-3,J)+2(j^+l)T*(N-2fJ)-(l+l/2(N-3))T*(N-ltJ)=D(N-Z) 
-(l-l/2(N-2))T*(N-2,J)+2(^+l)T*(N-l,J)  =  D(N-l) 
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where 

Dd)  =(^^^i,jH)^(f-c^)r(rJj;+(1^T(i,j+() 

(J=2,3 N-l)   (1=2, 3,..., N-l) 

and 

5(2)  =  D(2)  +  0.5;   D(N-l)  =  D(N-l)  +  1  +  l/£(N-2)). 
The  coefficients  needed  for  the  algorithm  subroutine,  TDIAG,  are 

B(I)  =  2(^+1)    (1=2, 3,..., N-l) 

A(I)  =  -(1-1/2(1-1))   (1*3,4 N-l) 

C(I)  =  -(1+1/2(1-1))   (I=2,3,...,N-2). 

Implicit  in  the  6  direction 

*[/ *^]T<-1^-  fsjjr^i  TCI, 3)  =  D(2) 

J^1(I'1^ t^Ll^mi)l^La,'L  _  =D(3) 

^J(zJ-i)^[f^^:]T(XJT)--^wlr{z1Tn)  =  D(J) 


l^Ta;^)^[^^(^Jr(r^-i)--L^T(i,-j-.)     =  d(n-2) 


where  D(J)  =   [l»^f]  T*fr/Jj)+4^)T*(tJj)+  [l  tj^j ]r*(l.H,;r) 

(1=2, 3,..., N-l)    ,      (J=2,3 N-l), 

and     D(2)   =  D(2)  +  (f^T^T  D(N-l)  =  D(N-l)    +     ^i  ^ 


with  coefficients 


B(J)  =  2{f+l/&I-lfK*))        (J=2,3 N-l),  (1=2,3 N-l) 

A(«D  -     TwFZF  U=3,^,..., N-l)  ,  (1=2,3 N-l) 

C(J)   =  A(J)  (J=2,3,...,N-2)  ,(1=2,3 N-l). 
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Note  that  in  these  equations,  each  coefficient  has  I  reduced 
by  one.   This  is  a  consequence  of  having  to  avoid  zero  subscripts. 

The  flow  charts  for  the  solution  program,  IAEMC,  and  the 
algorithm  subroutine,  TDIAG,  are  given  in  Figures  21  and  2. 

The  analytic  solution  is  given  in  Jaeger  Ll4j   f°r  ^e   points 
r  =  0.25,  0.5  and  0.75;  0  =  15°,  30°  and  kf  (the  temperature  is 
symmetric  about  the  ^5  degree  line)  at  times  between  0.01  and  0.2. 
These  temperatures  are  given  to  three  significant  figures ,  but  the 
author  states  that  calculations  were  made  to  four  significant 
figures  so  that  the  third  figure  is  considered  to  be  correct. 

Figure  22  shows  a  printout  of  the  intermediate  numerical  sol- 
ution, TST(I,J),  and  the  numerical  solution,  TEND(I,J),  for  time 
0.1.   For  this  solution  we  used  Ar  =  1/12,  A0  =  TT/24  and  At  =  0.005, 
Table  h   is  a  comparison  of  the  analytic  and  numerical  solutions  at 
time  0.1. 

Table  k.     Analytic  and  Numerical  Solutions  for  Time  0.1 


r— >                  0. 

25 

0.. 

5 

0. 

75 

Anal.    Sol. 

Num.    Sol. 

Anal.    Sol. 

Num.    Sol. 

Anal.   Sol. 

Num.    Sol. 

le 

15° 

0.979 

0.9788 

0.963 

0.9631 

0.957 

0.9574 

30° 

0.946 

0.9453 

0.906 

0.9053 

0.891 

0.8907 

*5° 

0.952 

0.9518 

0.917 

0.9166 

0.904 

0.9037 

Nine  different  solutions  were  obtained  using  combinations  of 
at,  a9  =  1/12,  TT/2^;  1/24,  TT/48;  I/36,  TV/72  and  At  =  0.005, 
0.0025,  0.00125  for  time  0.02.   Table  5  shows  these  results  at  point 
(0.5»  45  ).   The  difference  listed  in  this  table  is  the  difference 
between  the  numerical  solution  and  the  analytic  solution  at  this 
point.   The  analytic  solution  is  0.168, 
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Table  5.  Results  ->f  solutions  using  various  Ar.  A0,  and  •&  t 


at  the 

point   (0. 5. 

4<?°) 
(Ar)2+(  *0)2 

Num. 

At 

&r 

ao 

So  In. 

Difference 

0.00125 

0.0277 

0.0436 

0.00393 

0.167 

-0.001 

0.0025 

0. 0277 

0.0436 

0.00517 

0.168 

0 

0.005 

0.0277 

0.0436 

0.00767 

0.171 

0.003 

0.00125 

0.0416 

O.O654 

0.00728 

0.168 

0 

0.0025 

0.0416 

0.0654 

0.008 53 

0.169 

0.001 

0.005 

0.0416 

0.0654 

0.01103 

0.172 

0.004 

0.00125 

0.0833 

0.1309 

0.01345 

0.174 

0.006 

0.0025 

0.08^3 

0.1309 

0.01470 

0.175 

0.007 

0.005 

0.0833 

0.1309 

0.01720 

0.176 

0.008 

Figure  23  is  a  plot  of  the  absolute  difference  versus  At  + 

2        2 
(at)  +  (  A6)  using  the  data  given  in  Table  5.   The  line  shown  in 

this  figure  was  obtained  by  least  squares  and  has  the  equation 

0.626x  -  2.88  (we  have  scaled  all  values  by  a  factor  of  10  ,  for 

convenience).  Here,  as  in  Example  1,  the  data  appears  to  fit  the 

straight  line,  but,  as  previously  stated,  we  cannot  take  this  as 

verification  that  the  discretization  error  is  of  0  j_  At  +  (  •£  r) 

2 1 
+     (  A ©)  J  since  we  do  not  have  a  predicted  value  for  K.   The 

analytic  solution  as  given  in  [l4J   is 

u(r,0,t)  =  l-ff  £_         ,   Z_,  ^    —2 -ri  3 

where  J     is  the  Bessel  Function  of  order  s,    (s  =  2(2n  +  1)9, 
n  =  0,1,...),   and  0^         (w\  =  1,2,...)  are  the  positive  roots  of 

S  flu 

J  (CO  =  0.   Thus  we  see  that  it  would  be  very  difficult  to  obtain 

5 

the  derivatives  u. , ,  u    ,  and  u~,  .  needed  to  predict  a  value  for  K, 
tt   rrrr      0060         K 

The  complexity  of  the  above  equation  also  leads  us  to  conclude 
that  the  IAD  method  yields  a  simpler  program  than  one  involving  the 
analytic  solution. 
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Another  advantage  of  the  IAD  method  in  this  example  is  that, 
by  solving  the  equation  in  polar  coordinates,  we  avoid  the  necessity 
of  interpolation  at  the  circular  boundary.  Such  interpolation  be- 
comes necessary  when  placing  a  rectangular  net  on  a  circular  region, 
since  not  all  of  the  grid  points  will  fall  on  the  curved  boundary. 


59 


CONCLUSIONS 

From  the  numerical  results  of  the  three  examples,  we  have 
seen  that  accurate  results  are  possible  when  using  the  IAD  method. 
The  results  for  the  parabolic  problem  in  the  circular  region  suggest 
the  applicability  of  the  method  to  many  types  of  problems,  including 
fluid  flow  and  electric  potential.  When  the  differential  equation 
is  written  in  polar  coordinates, a  possible  difficulty  arises  when 
dealing  with  derivative -type  boundary  conditions  at  the  origin, 
since  the  differential  equation  has  a  singularity  at  that  point. 
(Temperatures  on  the  boundary  were  specified  in  Example  3t  thus 
avoiding  this  problem. )  Such  a  boundary  condition  would  require  that 
additional  approximations  be  made.   However,  the  method  seems 
ideally  suited  to  problems  involving  annular  regions.  As  pointed 
out  previously,  one  advantage  in  using  the  IAD  method  on  equations 
in  polar  form  is  that  there  is  no  need  for  interpolation  at  curved 
boundaries.  A  second  advantage  is  that  the  computer  program  for 
the  L.D  method  is  probably  simpler  to  write  than  a  program  to  obtain 
temperatures  at  the  grid  points  using  the  analytic  solution  (if  known), 

Since  Example  2  involved  an  elliptic  problem  over  a  rectangular 
region,  we  were  able  to  calculate  optimal  values  for  the  iteration 
parameter  V?  .   For  non-rectangular  regions  the  same  procedure  is 
not  applicable;  we  wdu*»,  therefore,  be  in  doubt  as  to  the  number  of 
iterations  required.   However,  it  is  possible,  when  using  simple 
values  of  y9  ,  to  gauge  the  progress  by  calculating  the  sum  of  the 
squares  of  the  errors  after  each  complete  iteration  to  see  what  the 
effect  of  a  particular  value  of  ^  might  be.   Details  for  this  pro- 
cess as  well  as  other  examples  of  solutions  in  non-rectangular 
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regions  may  be  found  in  Douglas  and  Peaceman  [ll]  .   In  addition, 
Wachpress  Ll6J  suggests  a  combination  of  the  IAD  method  and  an 
iteration  method  due  to  Lanczos  for  elliptic  problems  in  non-rect- 
angular regions.   Douglas   [15]  ,      [l7j   also  gives  details  for 
extension  of  the  IAD  method  to  mildly  nonlinear  problems  and  to 
problems  in  three  space  variables. 

Some  doubt  remains  whether  the  data  obtained  in  Examples  1 
and  3  actually  satisfies  the  predicted  error  bounds.  As  we  have 
seen,  the  data  obtained  from  the  numerical  solutions  is  not  sufficient 
in  itself  to  verify  predicted  results.  The  whole  question  of  dis- 
cretization error  is  found  to  be  much  more  complex  thar.  was  supposed, 
and,  in  fact,  could  be  the  subject  of  a  separate  paper. 

The  three  examples  given  show  that  the  IAD  method  involves 
essentially  the  same  equations  for  both  elliptic  and  parabolic 
problems.   The  algorithm  subroutine  used  to  solve  the  simultaneous 
systems  generated  by  the  IAD  equations  is  exactly  the  same  for  both 
problems  and  requires  only  a  change  of  input  parameters.   This  fact, 
and  the  satisfactory  results  obtained,  coupled  with  the  ease  of 
application,  leads  us  to  conclude  that  the  IAD  method  has  a  wide 
range  of  application  in  the  field  of  partial  differential  equations. 
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DIMENSION  A,B,C,D 
TST,  TEND,  TEMP,  OT, 
RHO,AN30L,  W,  2 


Ns  It 

NA=N-| 

XINC-O.I 


DO  S    J=  l,M 
DO  5    I  =  I,  M 


11 

TSTCr,J)=0 
TENDU,J)=0 

— 

— *- 

D0  6    JsJ.NA 
DO 6    I-  I.NA 

»- 

w 

TENDCIJ) 
r  l.O 

1». 

READ 
VALUES 
OF    DT 

IS]  RHO(n= 

XINC^XlNC 

DTCI) 

TAUSO 
L8    1 

TAUs 
TAU*DT(t) 

DO  15  1-1,12 

B(«)s2.0*RMO<K0 


F=2,0*RH0-».O 


DO  20  Xs  l,NA 


20J  A(X)s-I.O 
BCJ)-  B0> 

ccn=-i.o 


C(0  = 
-2.0 


DO  30  OH*  A 
DO 40  U  I,  MA 


fo 

0(1^2.0*  TEND  U, 2) 
t  f*  TEND  (I,») 

X  ' 

M)\ 

/*. 

Da)-TEND(I,J-h  + 
F-TENDa.JJ^TENDU.Ji'l) 

M 


CONTINUE 


CALL  TDIAfr 

CNAAft.C.D, 

TEMP) 


00  30X*I,NA 


30 


T5TU,J) 
sTEMPd) 


DO  50  Isi.MA 
D060  7el,MA 


-^ 


-»-<  ip«-i) 


D(J)s2.0«TSTa,J)  1^ 

+  F*TST(I,J) 


D(J)=TST(I-I,J)    * 
F#TSTa,J)  +  TST(I*l,J) 


FIGURE  1  FLOWCHART  FCR  PROGRAM  IADMT 


%> 


CONTINUE 

— ►- 

CALL  TDIAfr 
(NA,A,B,CiO,TEMP) 

0050  J-  i,NA 

50| 

TEND  CI, J) 

s  TEMPCJ) 

READ  VALUES 
or  AM  SOL 

DO  110  Isl,NA 

OOllO    Jr  |,NA 

Wtt,J>=ANSOlCr,J) 

-TEND  CI, J) 

TTol  ansolCI,J) 


PRINT 

Tir*ie,TAu 


fPRlNT  NUnffRlCAL? 
SOLUTION,  TEND 


fPRlNT  ANALYTIC  ^ 


SOLUTION,  ANSOL 


f         PRINT  \ 

DIFFERENCE,  W 


fPRiNT   PERCENT 
DIFFERENCE,  Z 


OTS  At 
XINC  e  A* 

RHOc     tf 


FIGURE   1      (COJTINUED) 
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SUBROUTINE  TD/AG- 
(NC,A,B,C,D,TEMP) 

DIMENSION    A.BjC.D, 
BETA,  GAMMA, TEMP 

BETA(Z^BU) 
GAMMA  U)  = 
D  CD /BETA  00 

BETA(Z)  «BCI) 

-au>*cci-o 

— *- 

D05  I=3,NC 

BETA(I-0 

CAMMAU)=LDU> 

-A(I)*frAMMACl-Ol 
-gj  /BETA  (I) 


TEMP(NC) 
=  GAMI»1A(NC) 

ND=  MC-Z 

• 

DO  10  K  =  I,ND 

IsMC-K 


TEMPCDs  6AMMACI) 
-CCI)«TEMP(I-H) 
BETA  CI) 


-HEND 


FIGURE  2  FLOWCHART  FOR  ALGORITHM  SUBROUTINE 
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DIMENSION       A,*,C,0 
TTND,T5T,  Tfc>IP,RHO, 

v,w,i 

Nell 
XINC  =O.I 
RfAO    V 

DO  2   Js  i,  N 

DO  X  1  =  1, M 

UTSTa,J)sO 
TENO<I,J)sO 

► 

srT 

VALUES 
OP  RMO(K) 

00  9  Jsl,N 
DO?  1=1, N 

YJs  J 
XI  =  I 

Ys(YJ-I.O)/iO.O 
Xs(XX-I.O)/iO.O 


£J  TST(I,J) 

x*Y/(i.o-y+x) 

TE-NOCr^rTSTCl.J) 


PRINT  TfeV4DU,jF 
Isl.N     Js  l,N 


Kff  I 


sr 


G0)sRHO*2,O 
FrRHO-X.0 


00  10  IeZ,N 


iOjA(I)r-I.O 
BCD  =6(0 
C(I)=-i.O 


A(f4)s  -2.0 

B(M)cRH042.0(i.O* 
XINC) 
NAsN-l 


DO  2.0   J*\tH* 
DO  90   IcZ.N 


20Jd<I)sTEND<I,J-i) 
♦  P^TrNOd,/) 
+  T£ND<X,T+0 


CALL    TDlAO 
C^iAjejC^D^TffMP) 

TST(IJ) 

=  TEMP(I) 

DO  20  t=2,N 

±o 


DO  40  IfZ,N 
00  SO  /=*,** 


^Fd-N) 


0U)s  i.O^TSKNA^) 
+  [«rtO-2.0*(«.OtXlNC)]*iTST(N,j) 


*0 


DM=TSTCm,J) 
tF*TST(I,T)  *TST(l>l,J) 


FIGURE  6      FLOWCHART   FOR  PROGRAM  IADM 
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M 


CONTINUE 


D<NA)  4-1,0 


call  tdiag- 
(naab^d.temp) 


40) 

TEN  OCX, D  * 
TEMPCT) 

DO  205"  J'Zfl4 
00  205    Is*,N 

0040  J=Z,NA 

->- 

vgalor)5Va,J)-nwp(r,j') 


"PRINT  RHO(K); 

ITERATION 

_  NUMgfRtK^ 


"PRINT  NUmSiCAL^ 
SOLUTION)  TEND 


PHI  NT   ANALYTIC 

Solution  ,  V 

v 


if*  INT  DlFFERCNC? 
W 


'PRINT   PERCENT  ^ 
DIFFERENCE,  t 

XINC  s   *% 
RMO  r  f 


FIGURE  6      (CONTINUED) 
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DIMENSION    A,6,C,D 
TSTtTEHD,T€MP 

N  r  13 
NA:  N-l 
RN  =  NA 

RlNC  s  i.O/ftrsJ 
ALFA  =  PJ/Z4,0 

AS  s  ALFA^ALFA 
OT^O.OO? 
TAUs  0.0 

RMO  3 
RlNC*RWC 
DT 

DOS  Jr!,N 
DO  f  1=  ',N 

— ►- 

Ls  1 

U  T5T<I,X)cO 
Tf*ID<I,,r)sO 

TffNOttjOvi.O 

DO  6  I* l,M 

fJ  TST(I,N)r|.o 
TCN0(X|N)8|.O 

TSTO»J>*  I.O 
TCNDOjttel.O 

DO  8  J=  '» N. 

8|  TSTfrjDci.O 
TENDC^DcLO 


ZJ  TAU  = 


TAU*DT 


OOIOJ;1,MA 


ajc  r-i 

6(7)- Z.0 

♦  (KHO-M.O) 

iSl  ccj)  = 

-(|.O*I.0/CZ.O*AJ)) 
A<J)c 

-(f.0-l.o/(x.0*A7)) 

DQZQfrZ,** 
DO  SO  I:  2,  NA 

AIS 

Ci-0«<r-i) 

DCI)s  0.0/(AI*A*))*T£ND(r,J-/) 
+  «,0«(RHO-I.O/<AZXAS»*TrNDCI,J) 

*<».0/CAX*A3»*TeN0(I,J>0 

FIGURE   21      FLOWCHART   FOR  PROGRAM   IADMC 
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CONTINUE 

OCX)  s  DUltO.f 
OfNA)  8     D(NA)+  1.0 
♦  I.O/(Z.O*(*M-/.0)) 

CALL  TDIAO 

CNA<A,slclo,reMO 

— ►- 

2oJ 
TSTCIJ) 

sTEMP(I) 

AIs(I-l)«(I-0 
All  sl-l 

00  20  I=Z,NA 

— ►- 

D040W,NA 

— *- 

^00  50  K5Z.NA 


12)    C(K)r-'.0/(AI*AS) 
B(K)-2.0*Cl.O/CAX*AS)  +  RHO) 

A(K)«  C<K) 


^00  60  J*2,NA 


OCJj  =  0.0-I.O/a.O*AIO)*T$TCl'M) 
+  2.O»CKMO-l,0)%TSTCX,«r) 

♦  0.0  +  l.O/r*.o*AX/))*TSTCX-H,.r) 

6£| 

continue: 

DU)  =  D(2>M.0/(Ar*AS) 
P(MA)=D<NAVh,0/W*AS) 

CALL    TDlAfr 
0tA,MjC#OjTY*») 

DO  4.0    T:7  WA 

[l£l 

TENDCI.J-) 
=  TEMPCT) 


/^"     PRINT 
TUMC^TAU 


/      PRINT 
INTERMEDIATE 
SOLUTION;  TST 


r         PRINT 

NUMERICAL 
SOLUTION,  TCHD 


RlNC  =  AV 
ALFAs  aG 
0T=  6* 
RHOr  f 


-^i 


LsL4  I 


FIGURE  21      (CONTINUED) 
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